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ABSTRACT 


A  sea,  swell  and  surf  program  is  improved,  tested  and  evaluated  on  a 
micro-computer  (HP-9845B).  Sea  swell  is  calculated  by  a  two  dimensional 
spectral  model.  The  energy  balance  equation  is  tested  for  different  cases  of 
wind  velocities  and  water  depths.  Satisfactory  agreement  is  observed 
between  the  offshore  model  and  expected  wave  heights  for  a  15  knot  wind, 
but  the  model  overbuilds  wave  energy  for  a  30  knot  wind.  Wave 

transformation  is  described  by  a  one  dimensional  random  wave  model  in  which 
the  wave  heights  are  described  using  the  Rayleigh  distribution.  The  obtained 
solution  of  the  random  wave  field  is  used  to  predict  the  longshore  currents. 
An  empirical  formula  for  determining  the  breaker  parameters  is  developed, 
based  on  beach  slope  and  incident  wave  steepness.  The  improved  model  is 
tested  using  an  undulated  bathymetry  to  validate  the  model  physics.  The 
model  outputs  of  wave  height  and  current  are  compared  with  data  acquired 
from  a  wave  tank  and  natural  beaches.  The  model  is  found  to  accurately 
forecast  wave  heights,  breaker  location,  breaker  type  and  longshore  currents 
for  several  sets  of  conditions.  Model  limitations  are  discussed  and 
recommendations  for  further  improvement  are  made.  fa  *  i  <ji)  f 'd  5  ;  ;  ~ 
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I.  INTRODUCTION 


A.  OPERATIONAL  REQUIREMENT  FOR  WAVE  FORECASTING 

Wave  and  surf  forecasting  has  been  an  important  part  of  Naval 
Amphibious  operations  since  the  beach  landings  of  World  War  II.  The  need 
for  the  prediction  of  ocean  and  surf  waves  for  military  purposes  led  to  the 
first  real  attempt  to  quantitatively  model  ocean  and  beach  waves.  Ocean 
engineers  and  physical  oceanographers,  such  as  M.P.  O'Brien,  H.U.  Sverdrup 
and  W.H.  Munk,  made  great  strides  in  sea,  swell  and  surf  modeling  during 
the  1940  -  1955  time  period  (Bascom,  1980).  The  wave  and  surf  predictions 
made  by  the  U.S.  Navy  are  still  based  on  these  works.  During  the 
intervening  time  period,  significant  advances  in  knowledge  of  ocean  wave 
dynamics  have  occurred,  and  it  is  time  to  update  prediction  techniques. 

Although  warfare  technology  has  advanced  along  all  fronts,  the  standard 
amphibious  beach  assault  is  still  a  viable  and  preferable  method  of  projecting 
forces  ashore  for  many  tactical  scenarios.  HQwever,  amphibious  assault  craft 
are  vulnerable  to  many  types  of  wave  dynamics.  In  general,  landings  are 
affected  by  the  initial  ocean  sea  and  swell  where  the  boats  are  launched  as 
well  as  by  the  breaking  waves  in  the  surf  zone,  along  with  the  accompanying 
longshore  currents.  Since  in  a  major  landing  tens  of  thousands  of  lives  are 
at  risk,  an  accurate  assessment  of  the  sea,  swell  and  surf  is  vital  to  a 
successful  amphibious  landing  operation  (Joint  Surf  Manual,  1976). 

In  addition  to  the  assault  landing  itself,  an  amphibious  operation  normally 
concludes  with  the  construction  of  a  temporary  harbor  facility  to 
accommodate  the  resupply  of  the  troops  ashore.  The  survivability  and  design 


characteristics  of  these  structures  depend  in  large  part  on  the  environment 


around  them.  The  ability  to  model  wave  height  and  direction,  and  longshore 
currents  is  necessary  to  operational  planners  and  coastal  engineers  in 
designing  and  emplacing  these  facilities. 

« 

In  the  past  twenty  years,  separate  models  have  been  developed  to  deal 
with  ocean  waves,  breaking  surf  and  longshore  currents.  In  1981,  the  U.S. 
Navy  contracted  to  have  a  model  developed  which  would  combine  three 
complementary  models.  The  goal  of  this  comprehensive  Sea,  Swell  and  Surf 
Program  (SSSP)  model  was  to  predict  the  wind  generated  wave  height, 
calculate  the  surf  zone  characteristics  based  on  the  wind  generated  energy 
propagation ,  and  use  a  two  dimensional  grid  to  compute  the  longshore  current 
field.  The  model,  developed  by  Wang  and  Chen  (1983),  was  written  in 
FORTRAN  and  designed  to  run  on  a  main  frame  computer.  To  facilitate  the 
model's  use  by  Naval  Oceanographers  in  the  fleet,  the  model  was  converted 
into  BASIC  and  implemented  on  the  HP-9845B/275  micro-computer 
(Devendorf,  1985).  The  primary  goal  of  this  thesis  is  to  modify  the  SSSP 
model  to  improve  its  computational  efficiency. 

B.  WAVE  THEORY  BACKGROUND 

The  water  surface  elevation  in  the  ocean  can  be  viewed  as  a  large 
number  of  individual  wavelets  at  various  frequencies.  The  waves  that  con¬ 
tribute  the  most  energy  to  the  ocean  wave  system  are  wind  generated  waves 
with  periods  of  1  to  30  seconds.  The  size  and  period  of  these  waves  are  a 
function  of  the  velocity  of  the  wind,  the  duration  of  the  storm  and  the 
distance  over  which  the  wind  blows  (fetch). 


V 
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The  simplest  theory  is  linear,  small  amplitude  wave  theory.  Sea  water 
is  assumed  to  be  incompressible  and  homogeneous,  and  surface  tension  forces 
are  assumed  to  be  negligible.  The  solutions  of  this  theory  describe  the 
behavior  of  the  individual  components  of  ocean  wind  waves  (Phillips,  1977). 
As  waves  enter  into  shallow  water,  wave  speed  varies  with  the  local  water 
depth  and  refraction  occurs. 

Initially,  refraction  of  monochromatic  waves  was  modeled  using  wave  ray 
theory  of  a  single  wave  (Arthur,  1949).  It  was  not  until  the  mid  1960's 
that  the  development  of  fast  computers  and  the  implementation  of  better 
physics  enabled  numerical  analysts  to  study  the  generation  of  wind  waves 
over  an  entire  frequency  spectrum.  Early  researchers  in  the  wave  spectral 
analysis  field  introduced  direction  variables  to  study  straight  and  parallel 
contour  wave  spectrum  transformation  (Karlsson,  1969)  and  included  the 
effects  of  bottom  dissipation  and  local  wind  generation. 

\s  computers  advanced,  more  complex  refraction  models  were  developed. 
\oda  et  at.  (1974),  usu  g  a  relaxation  finite  differencing  scheme,  solved  for 
the  stationary  wave  spectral  transformation.  The  model  that  is  under  study 
by  this  thesis  was  originally  developed  by  Wang  and  Yang  (1981)  and 
included  the  effects  of  bottom  friction.  Chen  and  Wang  (1983)  improved  the 
model  by  including  the  effects  of  non-stationary  waves. 

C.  CURRENT  MODELING  TECHNIQUES 

Two  numerical  modeling  techniques  are  currently  used  in  wave  and  surf 
forecasting  applications.  The  finite  element  method,  while  highly  complex, 
has  the  advantage  of  being  able  to  incorporate  irregular  localized  bathymetry 
and  has  stable  non-linear  convergence  properties.  Wu  and  Liu  (1985) 
introduced  the  method  in  their  study  of  wave  induced  nearshore  currents. 
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The  finite  difference  method  is  a  widely  used  numerical  method.  It  is 


based  on  a  linear  differencing  approach.  Finite  differencing  normally  uses  a 
rectangular  computational  grid.  This  method  was  used  by  Noda  et  al.  (1974  ) 
and  Shaiu  and  Wang  (1977)  in  their  studies  of  nearshore  circulation  and  wave 
energy  transformation,  respectively.  The  method  is  easier  to  implement  than 
the  finite  element  method  due  to  its  use  of  a  regular  grid  describing  the 
bathymetry,  and  is  employed  in  the  development  of  the  SSSP  (Wang  and 
Chen,  1983). 

D.  OBJECTIVES 

The  objective  of  this  thesis  is  to  improve  a  sea,  swell  and  surf  program 
which  is  implemented  on  the  HP  9845-B/275  micro-computer.  The  program  is 
divided  into  subroutines,  which  include  the  open  water  wave  generation 
module ,  wave  transformation ,  the  surf  zone  breaker  calculations  and  the 
longshore  current  routines.  To  test  the  model,  the  output  from  the  improved 
model  is  compared  with  data  bases  acquired  from  two  natural  beaches  as  well 
as  several  wave  tank  experiments.  A  users  manual  has  been  prepared  and  an 
attempt  has  been  made  to  interface  graphic  output  with  the  current  model. 
Limitations  of  the  model  are  discussed  and  future  improvements  suggested. 
This  model  is  proposed  for  use  by  operational  planners  and  Naval  Ocean¬ 
ographers  at  sea . 
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II.  OVERVIEW  OF  THE  SSSP  MODEL 


A.  OVERALL  PROGRAM  DESIGN 

The  Sea,  Swell  and  Surf  Program  (SSSP)  was  developed  by  Wang  and 
Chen  (1983).  The  SSSP  is  essentially  three  complementary  models  which  are 
merged  to  calculate  offshore  wave  height  fields,  surf  zone  information  and 
longshore  current  velocities.  The  three  modules,  ocean,  surf  and  longshore 
currents,  can  be  run  independently  or  the  offshore  wave  energy  can  be  used 
as  input  to  the  surf  and  current  models.  The  SSSP  is  a  numerical  model 
which  uses  finite  differencing  to  solve  the  governing  equations. 

The  original  version,  as  written  by  Wang  and  Chen,  was  written  in 
FORTRAN  and  designed  to  run  on  a  mainframe  computer.  The  version 
discussed  in  this  thesis  was  converted  to  BASIC  and  implemented  on  an  16 
bit  HP  9845-B/275  micro-computer  (Devendorf,  1985  ).  The  micro-computer 
version  consists  of  a  control  program  with  subroutines  to  calculate  wave 
direction,  height,  and  number,  and  output  surf,  swell  and  current 
information . 

The  SSSP  considers  an  offshore  wind  generated  energy  spectrum  which  is 
used  to  estimate  the  wave  growth  based  on  initial  inputs  of  wind  speed  and 
wave  height  or  spectral  energy  bands.  As  the  waves  traverse  down  the  user 
defined  grid,  changes  in  the  deep  water  energy  spectrum  due  to  bottom 
friction  and  irregular  bottom  changes,  are  taken  into  account.  If  the  surf 
zone  model  is  to  be  run ,  a  significant  wave  height ,  calculated  from  the  area 
under  the  wave  spectrum,  and  the  peak  frequency  are  used  to  provide  a 
monochromatic  wave  height  input  into  the  surf  model.  It  is  also  possible  to 
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input  monochromatic  wave  data  directly  into  the  surf  module  to  compute  the 


surf  zone  characteristics  without  running  the  computationally  intensive 
offshore  module  (Wang  ai.j  Chen,  1983). 


B.  SSSP  FLOWCHART 

The  SSSP  is  controlled  by  one  main  program  which,  in  turn,  can  call  any 
of  twelve  subroutines.  A  general  flowchart  of  the  entire  SSSP  shows  that 
the  SSSP  is  divided  into  two  main  modules  (Fig.  2.1).  The  user  is  prompted 
for  all  required  information  in  a  menu-driven,  interactive  fashion.  Initially, 
the  user  chooses  whether  to  run  the  offshore  or  the  surf  model.  General  and 
specific  input  is  requested  for  the  module  being  run.  The  open  water  or  surf 
condition  modules  are  calculated  and  the  proper  output  is  generated  by  either 
the  OUTSW  or  OUTSRF  routines. 

A  more  detailed  flowchart  of  the  open  water  condition  module  is 
presented  as  Figure  2.2.  OPEN  computes  the  wavenumber  as  a  function  of 
depth  and  then  loops  through  a  calculation  series  for  each  energy  band.  An 
adjustment  is  made  for  the  angle  of  the  swell/waves  with  respect  to  the  grid 
axis.  The  DIRECT  and  HEIGH2  subroutines  are  called  and  the  energy  from 
each  spectral  band  is  summed  and  converted  to  wave  height . 

The  original  SURFCON  subroutine  (Figure  2.3)  used  the  most  energetic 
energy  band  computed  by  the  OPEN  module  as  its  frequency  input,  or  used 
data  entered  directly  by  the  user.  The  wavenumber  was  computed,  the 
DIRECT  and  HEIGH2  subroutines  were  called,  the  average  wave  height  was 
calculated  and  a  series  of  breaker  height  and  locations  were  predicted. 

If  longshore  currents  were  desired,  the  NEAR  Cl  R  subroutine  was  invoked. 
The  model,  as  implemented  on  the  HP-9845,  used  an  iterative  process  on  a 
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Flowcharts  of  the  Open  water,  wave  Direction  and  wave 
Heiaht  subroutines  (OPEN,  DIRECT  and  HE1GH2). 
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two  dimensional  grid  to  solve  tor  the  mean  horizontal  velocities.  It  was 
found  that  this  scheme  was  computationally  unstable  as  well  as  requiring  as 
long  as  twenty-three  hours  to  converge  to  a  solution  (Devendorf,  1985). 

To  take  advantage  of  the  most  recent  advances  in  breaker  height 

modeling,  and  increase  the  computational  speed,  the  SURFCON  and  its 

ancillary  subroutines  were  revised  to  use  a  one  dimensional  random  wave 
model.  A' detailed  flowchart  of  the  new  SURFCON  module  (Fig.  2.4),  which 
replaces  the  original,  shows  the  simplifications  implemented  by  the  new 
model.  DIRECT2  is  a  simple  Snell's  law  calculation,  applicable  over  straight 
and  parallel  contours.  HEIGH3  uses  a  random  wave  model  which  incorporates 
the  probability  density  function  for  breaking  waves.  The  output  of  the  new 
SURFCON  module  includes  breaker  height,  number  of  breakers,  wave  direction 
and  surf  zone  width. 

LONSHOR  is  a  new  longshore  current  subroutine  (Figure  2.4),  which  was 
implemented  in  place  of  the  unstable  NEARCIR  to  simplify  and  speed  the 
calculation  of  the  longshore  currents.  This  model  uses  a  simple  radiation 

stress  balance  to  calculate  the  current  parallel  to  the  beach.  The  next 

section  will  discuss  the  physics  of  the  SSSP  modules  in  detail. 

C.  PROGRAM  MODULES 

1.  Offshore  Module 

When  the  SSSP  is  run,  a  menu  of  initial  choices  is  presented.  If  the 
user  chooses  the  "OPEN  WATER"  option,  the  OPEN  subroutine  is  invoked. 
This  offshore  module  computes  the  wave  height  and  direction  fields  as  a 
function  of  water  depth,  wind  speed  and  direction,  fetch,  and  the  initial 
energy  field.  The  initial  energy  can  be  described  by  a  single  significant 
wave  height,  direction  and  period,  or  as  an  implicitly  entered  energy 
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I NPUT2 

I NPUT 3 _ 

- - * - V 

y'N  EX  S  U  R  FCON  y 

COMPUTE  WAVE  NUMBER, 
CELERITY,  GROUP  VELOCITY 

CALL  01 RECT2 

CALL  HE  I GH  3 

CALC  THEORETICAL  MAX 
BREAKER  HEIGHT 

CALC  FIRST  BREAKER  LINE 
HE  I G  H  T  .  LOCAT  I  ON  AND 
DEPTH 

DETERMINE  BREAKER  TYPES, 
NUMBER  OF  BREAKER  LINES 


CALC  EFFECTIVE  SURF 


D  I  RECT2 


COMPUTE  WAVE  NUMBER  AND 
DEEP  WATER  CELERITY  C 

o 

DETERMINE  SNELL'S  CONST 
CALC  WAVE  ANGLES 
RETURN  SNELL'S  CONST 
RETURN 


J  H  E  1  6H  3 

COMPUTE  WAVE  NUMBER,  C 

AND  Cg 

CALC  INITIAL  ECg,  £  b 

CALC  INTERIOR  ECg,  £  b 

VALUES 

COMPUTE  RMS  WAVE  HEIGHTS 

CALL  LONSHOR 

RETURN 

* 

spectrum.  Water  depth  for  the  matrix  can  be  read  in  from  a  data  file  or 
entered  directly  from  the  keyboard.  A  tidal  height  correction  constant  can 
be  added  to  the  depth  field  to  take  into  account  the  state  of  the  tide  . 
The  bottom  type  can  be  specified,  which  determines  the  coefficient  of  friction 
for  the  wave  height  calculations, 
a.  Wave  Refraction 

Once  all  the  required  input  is  entered,  the  program  begins  its 
calculation  cycle.  The  first  field  that  is  calculated  is  the  wave  direction 
matrix,  using  a  steady  state  conservation  of  wavenumber  (Phillips,  2977): 


~(k  cos  6)  -  |^(k  sin  6)  »  0 


wnere : 


*  wavenumber 


2n  .  “1. 

"  T  (  } 


=  wavelength  (m) 

=  wave  direction  (clockwise  angle  from 
the  X  axis  to  the  wave  ray  (see 
Figure  7) ) 

=.  offshore  grid  spacing  (user  specified 
(see  Figure  5) ) 

=  longshore  grid  spacing  (user  specified 
(see  Figure  5) ) 


The  wavenumber  is  calculated,  for  the  frequency  band  being  considered,  using 


the  linear  dispersion  relation: 


(g  k  tanh  (k  h ) ) * ^ 


V."  *.*  ••  i  -V-V-V,'  V-.  ■.  •  v. 


where : 


c  =  angular  frequency  =  2r.f  (sec  ) 

-2 

g  =  acceleration  of  gravity  (m  sec  ) 
h  =  water  depth  (rn) 

f  =  frequency  =  ~  (sec  ) 

T  =  wave  period  (sec) 

Since  (2)  is  transcendental  in  k,  a  sixth  order  polynomial  fit  of  (2)  is  used 
(Hunt,  1979)  to  provide  a  faster  method  of  computing  k  than  the  iterative 
Newton's  method  used  in  the  original  Wang  and  Chen  (1983)  version  of  the 
SSSP.  The  wavenumber  k  is  given  as: 

k  =  - ^r-rrlS  +  (1  +  0.66667S  +  0.25550S2 

(gh) X/ 2 

+  0.16084S3  +  0.06320S4  -  0.021745* 

+  0.00654S6  +  0.00171S7  +  0.00029S8 
9  -1  1/2 

+  0.00011S  )  ]  (3) 


wnere : 


S  = 


2  ,h] 


Once  the  wavenumber  is  determined,  subroutine  DIRECT  is 
called  to  compute  the  wave  refraction  of  the  wave  components  over  the 
specified  grid.  DIRECT  starts  with  initial  values  of  the  wave  angles  for  each 
grid  point,  and  uses  a  finite  differencing  scheme  developed  by  Noda  (1972) 


to  solve  the  differential  equation  for  «  .  y  The  differential  equation  is 
center  differenced  in  the  offshore  ( X )  direction  and  forward  differenced  in 
the  longshore  (Y)  direction.  The  finite  difference  scheme  has  greater  weight 
on  the  forward  differencing  term  to  increase  computational  stability  (Abbott, 
1979  )  and  is  given  as: 

6,  .  =  sin  1{t~ — If  (k  sinfi.  ,  .  .  )  -  (1-2t)  (k  sin  9  .  ,  • ) 

±  '  J 


■IX 


)  -  cos6i|jtl: 


-  (*k  cos  6  ■  •  -i  )  ]  ) 

J-  /  J  ~  J- 


(4) 


where : 


t  =  weighting  factor  =  0.25 

A  first  guess  for  the  wave  angle  is  calculated  using  Snell's  law 

solution  for  straight  and  parallel  bottom  contours.  The  numeric  scheme  is 

iterated  until  r  .  .  is  within  0.005  of  the  previous  calculation.  The  10  X  13 
.  *  >J 

grid  point  system  used  by  the  SSSP  is  described  by  Figure  2.5.  The  initial 
differencing  scheme  used  by  the  model  is  a  simple  forward  stepping  scheme 
(Fig.  2.6).  The  convention  for  the  angle  -  is  illustrated  in  Figure  2.7. 
Theta  is  measured  with  respect  to  a  ray  drawn  perpendicular  to  the  beach. 
Therefore,  a  normally  incident  wave  would  have  a  Theta  value  of  180°. 
After  the  direction  matrix  has  been  specified,  the  subroutine  HEIGH2  is 


called  to  calculate  the  height  transformation  as  the  wave  moves  down  the 
grid  (towards  the  shoreline). 


Shoreline 


system  usnd  hy  SSSP 


2.  Wave  Height  Transformation 

This  subroutine  rakes  into  account  the  wave  direction  (previously 
calculated),  the  wind  generation,  and  the  bottom  dissipation.  HE1GH2 
employs  the  Phillips-Miles  growth  mechanism  for  the  wind  generation 
calculations.  An  initial  energy  field  is  developed  for  each  grid  point  by 
assuming  that  energy  is  conserved  over  straight  and  parallel  contours: 

—  (£  Cc  COS  r  )  =  0  > 5) 

where : 

E  =  wave  energy  (loules) 

Cg  =  group  velocity  (m  sec  1 ) 


The  initial  energy  field  is  derived  by  integrating  (5)  to  yield: 


Cc  cos  f 
~o  _ o 

o  Cg  cos  r 


(6) 


where  the  subscript  (o)  is  the  initial  condition  at  the  offshore  boundary. 

To  simplify  and  increase  the  model's  efficiency,  Chen  and  Wang  (1983  )  make 
four  assumptions: 

a.  Winds  are  steady 

b.  Wave  energy  within  a  spectral  band  is  restricted  to  that 
band  (Allows  superposition  of  wavelets). 
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c.  Each  frequency  component  can  be  described  by  a  single 
mean  direction ,  ~  . 

d.  No  wave-current  interactions. 

These  four  assumptions  allow  the  steady  state  energy  flux  for  each  spectral 
energy  band  to  be  written  as  ( Longuet-Higgins  and  Stewart,  1960;  1962  ): 

.  c  —  (  S  ;  f  )  Cg  cos  -  )  -  „c  —  IS  ;  f  i  Cc  sir,  ~ 

=  -■:-(  f  )  -  r  (  f  )  :  7  ) 

c  s 

where : 

S(f)  =  spectral  energy  flux 

=  density  of  fluid  (kgir.  ') 

tf)  =  energy  dissipation  due  to 

bettor,  friction 

•:  (f)  =  energy  generation  due  to 

wind  stress 

?  =  mean  wave  direction 

Tne  wind  generation  term  z  is  the  sum  of  the  linear  and  exponential 
contributions : 

e  =  c  _  c 
“m 


s 


(8) 


where : 


S_  =  Phillies  aeneration  =  cc; 
r* 


=  Miles  renerarron  =  pg0[S;f}] 


The  energy  generation  term  is  based  on  the  Phillips-Miles  growth  mechanisms. 
Starting  with  a  flat  ocean,  wave  growth  is  initially  linear  until  larger  waves 
are  present,  after  which  the  growth  is  exponential.  The  linear  growth  term 
(Phillips,  1957  )  is  based  on  random  atmospheric  pressure  variations  acting  on 
a  perturbed  sea  surface.  The  SSSP  first  calculates  the  pressure  fluctuation 
term  as  a  function  of  wavenumber  and  frequency: 


P(k,  c  )  = 


>-4„6 


*  (k  sin  : 


v2 


i  -  ( k  cos  :  -  ; 


rJ 


9) 


w  n  s  r  6  * 


W  =  wind  soeec  in  knots 


w 

1  o  p 

0 . 32.'.“ 


*  0.5  2/. 


,0.95 


=  ancle  between  wind  and  wave 


After  P(k,  z  )  is  calculated,  the  linear  growth  term ,  0  ,  is  determined  by: 


4n^kc  ^ 

a  =  -T-»  P(k,  c  ) 

"w? 


do: 
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where : 


t  =  density  of  water  (kg  m 

The  exponential  growth  term  becomes  important  after  waves  have  formed 
(LaBlond  and  Mysak,  1978  ).  When  flow  separation  occurs  in  the  lee  of  a 
wave  crest,  momentum  is  transferred  to  the  wave  because  of  the  pressure 
differential  resulting  in  exponential  wave  growth  (Miles,  1957;  1959a, b; 

1962) : 

£  =  (^~)  [  (|)  cos  6  -  b]  (11) 

where : 

a  =  constant  =  5.0 

S  =  ratio  of  air  to  water  density 

C  =  wave  celerity  (m  sec  ^) 

b  =  constant  =  0.90 

The  Phillips-Miles  wind  wave  generation  has  no  mechanism  for 
shutting  off  the  wave  growth.  In  nature,  wind  waves  will  grow  as  a 
function  of  wind  speed,  duration  and  fetch.  But  after  growing  to  a  certain 
size  and  steepness,  the  waves  will  become  unstable  and  break.  Therefore,  it 
is  necessary  to  have  a  limiting  function  so  that  the  waves  will  only  grow  to 
a  "fully  arisen"  sea.  Several  functions  have  been  proposed  to  introduce  an 
"equilibrium  value",  above  which  the  wave  growth  stops.  The  SSSP  uses  a 
Pierson  and  Moskowitz  (1964  )  fully  arisen  sea  in  deep  water.  The  high 
frequency  tail  is  described  by  a  Phillips  (1957)  equilibrium  condition, 
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modified  by  Kitaigorodskii  et  al.  (1975  )  and  Thornton  (1977)  for  intermediate 
or  shallow  water: 


where : 


S  ( f  ,  6 


B  g 2 f~5H(6) r (Wn) 


(12) 


H  ( 5 ) 


:(wn) 


w 


n 


8  4 

2  71  cos  i  accounts  for  angular 
spreading 

_-2 


2  w2  :  (w  )  -l 
(wj  [l  +  - — - - ] 


n 


sinh  [2W2Z  (Vf  )] 
n  n 


-(^)1/2 

“  V 


Z(W  )  meets  the  condition  that  Z  tanh(W  ~  )  =  1 
n  n 

■  9 

7  (W  )  =  1  in  deep  water  and  for  W  <  1 ,  7  ( W  )  — '  1/2  W  and  S  ( f ,  £  )  = 
n  n  n  n 

_  -»  _  9 

l/2qghdfJH(c)  where  q  is  on  the  order  of  10  "  (SSSP  uses  a  value  of 
0.073  ) . 


Bottom  dissipation,  c  ,  is  calculated  as  the  work  done  on  the 
bottom  due  to  the  bed  shear  stress.  Assuming  the  usual  quadratic  bed  shear 
stress  law,  the  wave  energy  dissipation  for  a  particular  spectral  frequency 
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band  is  described  by  Hasselmann  and  Collins  (1968): 


£d{f) 


pw  Cf 


U, 


cb(f) 


(13) 


where : 

=  bed  shear  stress  coefficient  =  0.01 

j ;  =  total  flow  field  bottom  velocity 

U^(f)  =  wave  induced  water  particle  bottom 

velocity  for  a  specific  spectral 
component 


Using  linear  wave  theory  transfer  functions  to  relate  the  surface  elevation 
energy  spectrum  to  the  bottom  velocity  field,  the  frequency  dependent 
dissipation  function  £^(f)  can  be  expressed  as: 


£d(f)  =  p  CfT(f)  S(f)  jUb|  (14) 

2  2  2  2 

where  the  transfer  function  T(f)  =  [g  k/o  cosh  (kh)].  The  total  wave 
induced  bottom  velocity  I |  is  calculated  by  integrating  the  spectral 
transfer  function  across  the  entire  frequency  spectrum: 

1/2 

:ub  =  I  /  T(f)  S(f)  df]  (15) 

0 

Bottom  dissipation  in  the  SSSP  is  only  calculated  when  deaLng  with 
intermediate  or  shallow  water  waves.  If  deep  water  wave  heights  are  being 
computed,  the  value  c  d  is  set  equal  to  zero  (Wang  and  Chen,  1983). 
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3.  Longshore  Currents 

One  of  the  most  serious  problems  with  the  original  model  was  the 
instability  of  the  nearshore  circulation  calculations.  As  much  as  23  hours  of 
CPU  time  were  required  for  the  subroutine  to  converge  to  an  answer.  It 
was  decided  that  a  one  dimensional  longshore  current  model  which  assumes 
straight  and  parallel  bottom  contours  would  be  substituted  for  the  unstable 
two-dimensional  model.  This  model  uses  a  simple  alongshore  radiation  stress 
balance  to  calculate  the  mean  longshore  current.  Although  no  circulation 
information  is  derived,  the  new  longshore  current  subroutine,  LONSHOR, 
provides  the  longshore  current  velocity  efficiently  and  accurately.  The 
physics  and  details  of  the  LONSHOR  subroutine  are  outlined  in  Chapter  II. 

4.  Surf  Zone 

The  original  SURFCON  module  made  calls  to  the  same  wave  direction 
(DIRECT)  and  wave  height  (HEIGH2)  subroutines  as  the  offshore  module. 
One  disadvantage  of  this  scheme  is  that  these  subroutines  are  very 
computationally  intensive,  especially  when  considering  the  narrow  width 
between  grid  points.  A  faster  wave  height  model  was  needed. 

The  state  of  the  art  in  breaker  height  models  uses  a  probability 
density  function  to  describe  breaking  wave  heights  (Thornton  and  Guza, 
1983).  Some  limitations  are  imposed  by  using  this  new  model,  however. 
The  model  assumes  that  the  beach  has  a  local  straight  and  parallel  bottom 
contour.  This  bathymetric  restriction  allows  longshore  features  such  as 
sandbars  to  be  modeled  accurately,  but  offshore  feature  such  as  channels  or 


canyons  are  smoothed  out  by  a  depth  averaging  subroutine.  It  is  felt  that 
this  assumption  is  reasonable ,  in  light  of  the  fact  that  many  beaches  of 


operational  interest  are  of  the  "straight  and  parallel"  type  within  the  first 
two  hundred  meters  from  shore. 

The  above  assumption  means  that  a  simple  Snell's  law  calculation  of 
the  wave  direction  is  possible.  A  short  subroutine,  called  DIRECT 2,  is 
implemented  to  calculate  the  wave  direction  matrix.  The  use  of  the  new 
height  and  direction  modules,  HEIGH3  and  DIRECT2  have  increased  the 
efficiency  of  the  SURFCON  calculations  immensely.  The  average  run  time  for 
the  surf  zone  model  predicted  wave  heights  has  decreased  from  over  four 
minutes  to  less  than  30  seconds. 
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HI.  MODEL  IMPROVEMENTS 


A.  OFFSHORE  WAVE  HEIGHT  GROWTH 

The  current  version  of  the  SSSP,  converted  to  run  in  HP  BASIC  on  the 
HP-9845B/275  micro- computer,  gives  reasonable  values  for  wave  height 
growth  for  wind  speeds  of  15  knots  or  less  (Devendorf,  1985).  However,  for 
higher  wind  speeds  of  30  knots,  the  SSSP  exhibits  a  growth  rate  that  is 
significantly  faster  than  the  growth  of  the  JONSWAP  spectrum  (Bishop, 
1983),  used  for  comparison  in  this  analysis.  The  JONSWAP  spectrum  was 
chosen  as  a  basis  for  comparison  because  the  JONSWAP  curves  were 
experimentally  derived  from  the  limited  fetch  North  Sea  area.  It  is  foreseen 
that  the  SSSP  model  may  be  applied  to  similar  operational  areas. 

For  30  knot  winds,  the  model's  wave  growth  not  only  occurs  too 
quickly,  reaching  a  maximum  value  within  50  nautical  miles,  but  the 
magnitude  of  the  growth  is  too  high  (maximum  wave  height  of  6.0  meters 
compared  to  a  JONSWAP  maximum  of  4.2  meters.  Before  sensitivity  tests 
were  run  on  the  wind  generation  portion  of  the  model,  the  code  was  checked 
carefully  to  ensure  that  the  equations  were  being  implemented  correctly. 

To  determine  the  sensitivity  of  the  model  wave  growth,  several  runs 
were  made  with  changes  to  the  key  growth  parameters.  These  parameters 
include  A,  the  Phillips  linear  growth  term.  Equation  9  (Phillips,  1957);  B, 
the  Miles  exponential  growth  term,  Equation  10  (Miles,  1957;  1959a, b; 

1962),  and  R,  a  weighting  factor  set  equal  to  1.0  in  the  original  model. 


7.  r.  -r.  «rv* 


The  first  run  was  used  as  a  standard  to  compare  the  modifications 
against.  The  inputs  for  the  OPEN  module  were: 

Grid  Spacing  X  =  Y  =  10  NM 

Wave  Direction  Theta  =  180°  (normal  incidence) 

Wind  (W)  =  30  and  15  kts  ®  270° 

Shoreline  =  270° 

Bandwidth  (Df)  =0.01  Hz 
Frequency  bands  =  0.05  -  0.14  Hz 

The  output  for  these  "standard'1  model  runs  are  shown  as  Figures  3.1  and 
3.2.  The  overbuilding  of  the  open  water  waves  for  W  =  30  kts  is  evident 
when  compared  with  the  JONSWAP  spectrum  wave  growth  (Fig.  3.1).  For  a 
more  moderate  case  where  W  =  15  kts,  the  model  wave  growth  shows 
satisfactory  agreement  with  the  JONSWAP  growth  (Figure  3.2).  To  see  the 
impact  of  changing  the  weighting  factor,  R,  two  runs  were  made  setting  R 
equal  to  0.5  and  0.1  respectively.  The  predicted  wave  height  with  R  equal 
to  1,  0.5  and  0.1.  is  shown  In  Figure  3.3.  Setting  R  equal  to  0.1  (chain 
dashed  line)  gives  a  good  fit  for  high  wind  speeds  but  underbuilds  the  wave 
field  for  lighter  winds  (Fig.  3.4).  Clearly,  a  simple  change  of  the 
weighting  factor  R  will  not  resolve  the  problem. 

The  Miles  exponential  growth  term,  B,  is  in  the  denominator  of  the 
energy  calculation  so  its  value  was  increased  to  determine  the  growth 
retarding  effect  it  would  have.  Two  runs  were  made  with  B  doubled  (chain 
dotted  line)  and  with  B  increased  by  a  factor  of  10  (not  shown),  depicted 
by  Figure  3.5  The  result  of  increasing  the  exponential  growth  term  is  to 
flatten  out  the  upper  portion  of  the  growth  curve  which  is  not  the  correction 
that  needs  to  be  applied. 
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Model  output  (solid)  versus  JONSWAP  curve  (dashed) 
=  30  kts.  The  SSSP  substantially  over  develops 


Model  output  (solid)  versus  JONSWAP  curve  (dashed)  for 
=  15  kts.  For  low  wind  speeds,  there  is  good  agreement 
between  JONSWAP  and  SSSP. 


gure  3.4.  Sensitivity  test;  R  =  0.1  (chain-dot)  for  W  =  15  kts. 
Since  each  is  multiplied  by  R,  a  reduced  R  value  will  decrease 
the  low  wind  speed  case  as  well  as  the  high.  The  unchanged 
SSSP  model  (solid)  and  JONSIVAP  (dashed)  are  also  shown. 


Because  the  model's  rapid  growth  appears  to  occur  in  the  linear  region 
of  the  growth  curve,  the  Phillips  linear  term,  A,  was  analyzed  closely  to 
see  if  a  simple  constant  multiple  would  bring  the  wave  growth  in  line  with 
the  JONS  WAP  curves.  A  constant  multiple  of  0.125  is  applied  to  the  A  term 
(Fig.  3.6).  This  method  gives  a  better  fit  to  the  JONSWAP  curve  for  high 
wind  speeds,  although  the  SSSP  still  over  builds  the  wave  heights.  To  test 
this  multiple  at  lower  wind  speeds,  a  model  run  was  made  for  a  wind  speed 
of  15  knots,  for  a  bandwidth  of  0.05  Hz  and  frequencies  of  0.05,  0.10, 
0.15,  and  0.20.  The  results  of  the  model  run  are  compared  with  the 
equivalent  JONSWAP  curve  and  the  curve  from  the  existing  model  (Fig.  3.7). 
The  fit,  while  not  exact,  is  reasonable,  suggesting  that  the  A  term  is 
responsible  for  the  overbuilding  of  the  model  wave  heights. 

The  actual  cause  of  the  overbuilding  of  wind  waves  by  the  SSSP  is 
unknown.  Several  ad  hoc  "fixes"  to  the  model  are  suggested  with  the 
recommendation  that  further  research  on  this  problem  be  conducted.  A 
simple  constant  multiple  of  the  linear  growth  term,  on  the  order  of  0.10  to 
0.15  is  one  recommended  empirically  derived  modification  to  the  wind 
generation  calculation.  Another  possible  modification  is  to  have  the 
weighting  function,  R,  decrease  as  the  wind  increases,  to  act  as  a  damper 
on  the  over  generation.  At  present,  the  SSSP  wave  generation  routine 
remains  in  its  original  form.  It  Is  felt  that  the  model  is  adequate  for  low 
wind  speed  generation.  Operational  planners  should  be  cautioned  about  the 
high  wind  speed  generation  problems. 

B.  WAVE  REFRACTION  MODULE 

To  test  the  SSSP  model's  stability  during  the  wave  height  prediction  in 
the  coastal  zone,  a  sinusoidally  varying  bottom  was  explicitly  entered  using 
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Figure  3.6.  Sensitivity  test;  A  term  reduced  by  constant  multiple  of 
0.125  (clmin-dot)  tor  W  =  30  kts.  The  unchanged  SSSP 
model  (solid)  and  JONSWAP  (dashed)  are  also  shown. 


Figure  3.7.  Sensitivity  test;  A  =  0.125  A  (chain-dot)  for  W  =  15  kts. 

The  unchanged  SSSP  model  (solid)  and  JONSWAP  (dashed) 
are  also  shown . 


the  following  equation: 


h,  ,  =  0.025X  +  ch.  (r^-)  (1  -  (tt-))cos(^)  (16) 

(x,y)  d  Xb  Xfa  Liy 

where : 

z  =  bottom  perturbation  term  (0.2) 

hb  =  breaker  depth  (2.5  m) 

X,  =  offshore  distance  to  breaker  line 
°  (100  m) 

1.^  =  wave  length  of  bottom  perturbation  (  240  m) 

This  equation  provides  a  smoothly  varying  bottom  profile  in  both  the  X 
(offshore)  and  Y  (longshore)  directions.  To  check  that  the  bottom  profile 
was  being  accepted  correctly  by  the  model,  and  that  the  calculated  results 
were  reasonable,  approximately  25  model  runs  were  completed  with  varying 
input  parameters.  Selected  cases  are  presented  herein  for  discussion.  In 
the  cases  presented,  the  offshore  waves  are  normally  incident  to  the 
shoreline.  Obliquely  incident  model  runs  were  made  with  Theta  ranging  from 
+10  to  -10  degrees  from  normal,  for  comparison  with  the  normally  incident 
case . 

1*  Model  Experiments 

The  input  parameters  in  the  model  are:  significant  wave  height,  Hg 
(1  m);  wave  period,  T  (12  s);  offshore  boundary  depth,  h(l2-  m);  artificial 
viscosity,  T  (0.0);  and  bottom  perturbation  term,  £  (0.2).  The  beach 
faces  270°  and  the  normally  incident  wave  rays  come  from  270°.  This  gives 
a  Theta  value  of  180°  according  to  the  model  convention  shown  in  Figure 


The  model  results  for  wave  heights  and  directions  are  shown  in 
Figures  3.8-3.15.  The  longshore  depth  profiles  for  locations  inside  (X=3) 
and  outside  (X=8)  the  surf  zone  are  chosen  to  examine  the  wave  refraction. 
Grid  point  X  =  3  implies  a  distance  of  40  meters  offshore,  well  inside  the 

breaker  line,  while  X  =  8  is  140  meters  offshore  and  outside  the  breaker 

line.  The  depth  profile  presented  in  Figure  3.8  shows  that  the  depth  varies 
sinusoidally  in  the  longshore  direction.  The  beach  profile  along  grid  point  Y 
=  3  is  shown  as  Figure  3.9.  As  shown,  the  depths  are  slightly  deeper  than 
the  plane  beach  inside  100  meters  (breaker  position),  and  shallower  outside 
100  meters. 

As  waves  enter  the  surf  zone,  their  increase  in  height  is  a  function 
of  water  depth.  For  the  points  outside  the  breakers,  the  largest  waves  are 
expected  to  occur  where  the  water  depth  is  the  shallowest  due  to  shoaling. 
As  can  be  seen  from  Figure  3.10,  the  peak  waves  (dashed)  are  towards  the 

edges  of  the  domain,  but  they  do  not  coincide  with  the  minimum  in  depths  at 

0  and  240  meters.  This  anomaly  is  a  result  of  the  boundary  conditions  of 
the  model.  Also  note  the  asymmetric  peaks  at  40  and  200  meters.  This 
implies  errors  in  wave  height  prediction  due ,  in  part ,  to  the  wave  direction 
calculations.  These  errors  will  be  addressed  later  in  this  discussion. 

Wave  direction  is  calculated  in  terms  of  Theta.  An  angle  greater 

than  180  tends  to  turn  the  wave  ray  to  the  left  while  Theta  less  than  180 

turns  the  ray  to  the  right  (See  arrows  in  Figure  3.11).  The  model  at  least 

quantitatively  predicts  the  wave  refraction  pattern  correctly  for  all  cases, 

consistently  turning  the  rays  to  "attack"  the  shallower  water  depths  outside 

the  breaker  line.  The  model's  boundary  conditions  force  0^  =  9^  and  G  n_^ 

=  9  where  S  is  the  wave  direction  at  the  n**1  long  shore  grid  point.  This 
n 
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Figure  3.8.  Longshore  depth  profile  for  X  =  3  (solid)  and  X 
(dashed).  The  sinusoidally  varying  bottom  profile 
used  as  a  depth  field  for  testing  the  SSSP. 


a  w 
o 
2 
o  < 
■o  E-« 
^  CO 


.o  M 

®  « 


o  CO 
©  bn 


o*o  ro-  Z'O-  e*o-  ro- 

(j\t)  iviiKaaaiiia  Hxaaa  3aym 


45 


Figure  3.9.  Offshore  depth  profile  for  Y  =  3  (solid).  Depths  are 
deeper/  shallower  than  a  plane  beach  profile  of  slope  0.025 
(dashed) . 


Figure  3.10.  Wave  Heights  at  X  =  3  (solid)  and  X  =  8  (dashed).  Note 
asymmetric  peaks  at  40  and  200  meters.  Wave  heights 
inside  the  breaker  line  are  related  to  the  bottom  depth 
only . 
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accounts  for  the  departures  of  the  wave  direction  at  the  edges  of  the 
domain.  The  dotted  lines  are  the  interpretation  of  how  the  wave  direction 
patterns  should  behave.  In  the  middle  of  the  domain,  where  both  cases 
exhibit  a  nodal  point,  the  wave  direction  passes  through  180  or  normal  to 
the  beach  as  expected. 

The  final  test  conducted  compared  the  effect  of  Tau  on  the  wave 
direction  calculations.  Tau  is  an  artificial  viscosity  parameter  which  is 

intended  to  be  a  smoothing  term  for  the  calculation  of  wave  angle. 
Changing  Tau  from  0.00  to  0.25  had  a  negligible  effect  on  the  wave  height 
calculations,  but  the  difference  in  direction  increased  with  decreasing 
distance  from  shore.  The  case  where  X  =  8  is  plotted  in  Figure  3.12.  The 
introduction  of  the  Tau  term  (dashed)  decreases  the  magnitude  of  the  wave 
angle,  but  the  maximum  (minimum)  values  are  not  exactly  in  phase  with  the 
Tau  =  0.00  case.  The  maximums  are  shifted  slightly  toward  the  edges  of  the 
domain.  This  is  the  expected  behavior  of  the  smoothing  term. 

The  dissipative  model  exhibits  unexplained  behavior  inside  the 

breaker  line  as  depicted  by  Figures  3.13  -  3.15.  At  X  =  5  (Fig.  3.13), 
the  Tau  model  reduces  the  magnitude  of  the  wave  direction  field.  At  X  =  4 
(Fig.  3.14),  spurious  asymmetry  appears  in  the  direction  calculations  with 
the  Tau  term  applied.  Furthermore,  at  X  =  3  (Fig.  ■»  li>),  the  effect  of  Tau 
is  to  increase  the  magnitude  of  the  wave  direction.  This  apparent  angular 
deviation  has  no  effect  on  the  model's  wave  height  calculations  because, 

inside  the  breaker  line,  wave  height  becomes  a  function  of  depth  only. 
However,  this  instability  does  effect  the  longshore  current  calculations, 

which  are  sensitive  to  local  wave  direction  inside  the  surf  zone. 
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Wave  Directions  at  X  =  5  for  Tau  =  0.0  (solid)  and  0.25 
(dashed).  The  expected  modifying  behavior  of  Tau  is 
exhibited . 


Wave  Directions  at  X  =  4  for  Tau  =  0.0  (solid)  and  0.25 
(dashed).  As  the  wave  approaches  the  shoreline,  the 
Tau  modification  shows  irregular  behavior. 


'igure  3.15.  Wave  Directions  at  X  =  3  for  Tau  =  0.0  (solid)  and  0.25 
(dashed).  Contrary  to  its  intent,  close  to  the  shoreline, 
Tau  actually  increases  the  magnitude  of  the  wave 
direction. 


Another  indication  of  a  numeric  error  in  the  wave  direction 
calculation  was  noted  in  the  exact  center  (Y  =  7)  of  the  domain.  Since  the 
depth  profiles  exhibit  a  maximum  (minimum)  value  at  the  domain  center, 
theory  suggests  that  the  waves  remain  normally  incident  ••  t'tese  nodal 
points.  Therefore,  the  Theta  value  should  be  180  degrees  at  the  center  of 
the  domain.  The  model  output,  both  with  and  without  Tau,  exhibits  a 
monotonically  increasing  error  as  seen  in  Figure  3.16. 

2.  Improvement  of  the  Numerical  Scheme 

In  an  effort  to  eliminate  the  numerical  instability  in  the  wave 
direction  module,  the  finite  differencing  scheme  was  analyzed  to  see  if  an 
improvement  could  be  made.  The  original  scheme,  pictured  in  Figure  2.6, 
uses  a  simple  centered  and  forward  differencing  method  and  is  fairly  robust  in 
the  longshore  ( Y )  direction  but  is  less  so  in  the  offshore  ( X )  direction .  To 
increase  the  accuracy  of  the  numerical  scheme,  the  method  pictured  in 
Figure  3.17  was  substituted  in  the  DIRECT  subroutine  for  the  original 
scheme.  This  particular  method  is  a  well  accepted  technique  which  has 
proven  to  be  stable  and  free  of  damping  (Haltiner  and  Williams,  1980). 

Several  model  runs  were  made  to  test  the  new  technique.  Using  the 
standard  input  parameters  and  setting  Tau  equal  to  0,  the  first  result  was  a 
gain  in  accuracy  of  the  wave  direction  calculations.  The  domain's  central 
directions  were  exactly  180°  (to  three  significant  figures),  as  seen  in  Figure 
3.18,  and  the  refraction  angles  were  symmetric  about  the  center  of  the 
domain  as  expected  due  to  the  symmetrically  varying  bathymetry  (Figure 
3.19). 

A  second  result  of  the  implementation  of  the  improved  numerical 
technique  is  the  stabilizing  effect  on  the  wave  height  calculations.  The 
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OLD  SCHEME 


Figure  3.17.  A  more  robust  differencing  scheme  was  substituted  for  the 
original  numeric  scheme  seen  in  Figure  2.6. 


Figure  3.18.  Increased  stability  in  Theta  calculations  resulting  from  the 
new  numeric  scheme  (comgare  with  Figure  24).  The 
expected  value  of  180.000  is  achieved  with  the  new 
scheme . 


Figure  3.19.  Wave  Directions  for  Tau  =  0.0  at  X  =8  (solid)  using  tlte 
new  numeric  scheme  compared  with  the  old  scheme 
(dashed).  Dotted  lines  show  expected  behavior  at  the 


previous  asymmetric  wave  heights  are  now  symmetric  with  respect  to  height 
and  peak  location  (Fig.  3.20). 


3.  Boundary  Conditions 

The  SSSP  model  imposes  no  flux  boundary  conditions  on  the 
refraction  and  wave  height  calculations  so  that  the  left  and  right  boundaries 
are  specified  as  follows: 

Q(M,1)  =  Q(M,2) 

Q(M,N)  =  Q(M,N-1) 

where:  Q  =  the  quantity  being  calculated 
M  =  offshore  row  number  (10  to  1) 

N  =  the  last  longshore  column  (13) 

The  offshore  row  is  completely  specified  by  the  initial  conditions.  The 
calculations  then  proceed  from  offshore  to  onshore  and  from  left  to  right  on 
the  grid. 

No  flux  boundary  conditions  account  for  the  unusual  behavior  of  the 
wave  height  and  direction  curves  at  the  edges  of  the  domain.  To  satisfy  the 
condition  of  well-posedness  for  the  applied  model,  cyclic  boundary  conditions 
are  substituted  such  that: 

Q(M,l)  =  Q(M,N-1)  (17) 

Q(M,N)  =  Q(M,2) 

A  plot  of  wave  height  versus  longshore  distance  at  X  =  8  using  the  same 
inputs  as  the  previous  case  is  presented  in  Figure  3.21.  Note  the  shift  of 
the  peaks  toward  the  edges  of  the  domain  as  well  as  the  abrupt  increase  in 
the  magnitude  of  the  wave  height.  These  numerical  results  now  match  quite 
well  with  theory  which  requires  that  the  wave  rays  turn  toward  the 


shallowest  water,  implying  that  the  largest  wave  heights  should  occur  at  the 
edges  of  the  domain. 


The  wave  direction  calculations  also  show  improvement.  Close 
agreement  between  the  cyclic  boundary  condition  values  and  the  expected 
values  (dotted)  is  shown  in  Figure  3.22.  Additionally,  the  average  value  on 
the  boundary  segment  agrees  closely  with  the  exact  value  (solid  circles). 

4.  Results 

Based  on  the  results  of  many  model  experiments,  the  new  model 
appears  to  offer  an  increase  in  stability  and  computational  accuracy.  The 
improved  wave  direction  subroutine  consistently  turns  the  wave  rays  in  the 
correct  direction.  The  magnitude  of  the  direction  changes  is  qualitatively 
correct.  The  wave  height  subroutine  works  well  both  inside  and  outside  the 
surf  zone.  Wave  heights  are  highest  at  the  shoal  and  lowest  at  the  trench. 

One  drawback  to  the  original  wave  height  subroutine  is  the 
computational  speed.  When  running  the  program  to  generate  an  energy 
spectrum  for  sea  and  swell,  the  model  takes  approximately  45  minutes  per 
frequency  band.  Even  the  monochromatic  surf  module  is  slow  (4.5  -  5 
minutes)  because  of  the  number  of  iterations  required  to  converge  to  the 
preassigned  criterion. 

Least  acceptable  of  all  is  the  Nearshore  Circulation  module  which 
proved  to  be  unstable  as  well  as  requiring  as  much  as  23  hours  to 
asymptotically  approach  an  answer.  To  provide  fast,  accurate  longshore 
current  information,  a  one  dimensional  model  is  substituted  for  the  two 
dimensional  circulation  model.  While  some  potential  rip  current  information 
will  not  be  provided  by  this  simple  model,  the  one  dimensional  model  has  the 


Figure  3.20.  Wave  Heights  at  X  =  8  using  new  numeric  scheme  (solid). 

Note  expected  wave  height  symmetry  compared  to  original 
heights  (dashed). 


Figure  3.21.  Wave  Heights  at  X  =  8  with  cyclic  boundary  conditions 
imposed  (solid)  compared  with  the  wave  height  field  using  the 
old  boundary  conditions  (dashed)  The  wave  peaks  have  shifted 
to  the  edges  of  the  domain  as  suggested  by  theory. 


advantage  of  being  fast,  accurate  and  able  to  provide  the  most  important 
input,  the  longshore  current  distribution. 

A  state  of  the  art,  one-dimensional  wave  height  model  is  used  as 
input  to  the  longshore  current  model.  The  new  wave  height  and  longshore 
current  subroutines  are  described  below. 

C.  ONE  DIMENSIONAL  WAVE  HEIGHT  CALCULATION 

To  Increase  the  SSSP's  computational  efficiency  while  maintaining 
reasonable  accuracy,  a  one  dimensional  random  wave  height  model  was 
substituted  for  the  original  two  dimensional  monochromatic  model.  This 
model  characterizes  the  transformation  of  the  wave  height  probability  density 
function  (pdf)  from  offshore  to  the  shoreline  (Thornton  and  Guza,  1983). 
The  model  assumes  that  the  waves  are  narrow  banded  in  frequency  and 
direction  such  that  the  waves  are  specified  by  a  single  mean  frequency  f  and 
mean  direction  9  ,  that  the  bottom  contours  are  straight  and  parallel,  and- 
that  the  wave  conditions  are  stationary.  This  probabilistic  approach  to  the 
wave  breaking  problem  was  originally  suggested  by  Collins  (1970)  and  Battjes 
(1972).  The  form  that  is  used  here  is  based  on  integrating  the  energy  flux 
balance  equation.  The  calculated  wave  heights  are  a  function  of  the  shoaling 
conditions  of  the  integral  path  from  offshore  to  the  shoreline  ( Battjes  and 
Janssen,  1978). 

For  straight  and  parallel  contours,  the  energy  flux  is  balanced  by  wave 
breaking  and  the  frictional  dissipation: 


3  (E  Cq  cos  0 ) 


(18) 


where : 


<efa>  =  dissipation  due  to  wave  breaking 

<c^>  =  frictional  dissipation 

In  the  domain  of  interest,  it  has  been  shown  that  the  frictional  dissipation  is 
less  than  3%  of  the  dissipation  due  to  breaking  waves  (Thornton  and  Guza, 
1983)  and  will  be  neglected. 

The  waves  in  the  new  model  are  described  by  the  Rayleigh  wave  height 
distribution : 


P(H) 


2H  _  rms 


(19) 


where: 


H _ _  *  root  mean  scruare  wave  height 

where  Hrms  is  the  root  mean  square  wave  height.  This  distribution  takes 
into  account  the  fact  that  there  is  a  distribution  of  wave  heights  throughout 
the  entire  field  for  a  given  frequency.  Previous  random  wave  models 
truncated  the  Rayleigh  distribution  in  various  ways.  Recent  observations 
have  found  that  even  after  the  waves  break,  the  Rayleigh  distribution  is  still 
valid  (Thornton  and  Guza,  1983).  Therefore,  in  the  improved  model,  the 
waves  are  everywhere  Rayleigh  distributed.  The  probability  of  a  wave 
breaking  p^fH)  is  described  by  a  simple  weighting  of  the  Rayleigh 
distribution : 

Pb(H)  =  W(H)  p(H) 
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The  weighting  function,  W(H),  which  makes  larger  waves  at  a  given  depth 
more  likely  to  break,  is  given  by  (Thornton  and  Guza,  1983): 


H  4 

W(H)  =  l~jp] 


(20) 


wnere : 


7  =  adjustable  parameter  based  on 

field  data 

The  model  assumes  that  breaking  waves  are  similar  to  periodic  bores 
(LeMehaute',  1962)  as  depicted  in  Figure  3.23.  The  bore  dissipation  per 
unit  area  is  calculated; 


cbore 


(h2-h1) 

nin2 


3 

-0 


(21) 


where : 

Q  =  volume  discharge  across  bore 

B  =  breaker  coefficient  (described 
from  data) 


The  volume  discharge  parameter  is  described  for  a  linear  periodic  bore 
(Hwang  and  Divoky,  1970)  by: 

Q  =  Ch/L  (22) 

where  C  is  the  wave  celerity.  Using  linear  wave  theory,  the  wave  energy 
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and  group  velocities  are  given  by: 


E  =  pg  J  H2  p  (H)  dH 
0 


H 


rms 


(23) 


and, 


C 

gx 


2  k  h 


sirin  2kn 


cos  r 


(24) 


After  integration  and  the  substitution  of  the  probability  density  function , 
Equation  21  leads  to: 


where: 

i  =  mean  frequency 

This  equation  for  the  bore  dissipation  can  be  tuned  to  the  data  by  varying 
the  breaker  coefficient  B ,  ai.  .  the  adjustable  parameter  >  .  Substituting  the 
bore  dissipation  function  (25),  the  energy  flux  equation  (23)  is  given  by: 


■A 


^  g 


fB3  7 

4,5  rms 

y  h 


(25) 


d  .1  n  u2 
dX ( 8  p  9  H 


C  cos  6) 
rms  g 


Te  V  T  09 


fB3 


H 


(26) 


rms 


For  shallow  water,  the  group  velocity  Cg  is  approximately  equal  to  the  wave 

1/2 

celerity  C  which  equals  (gh)  .  Equation  (26)  was  chosen  to  be  used  in 
the  wave  height  module  because,  although  a  more  complicated  expression 


gave  a  slightly  better  fit  to  the  data,  Equation  (26)  provides  an  analytical 


solution  which  allows  a  detailed  study  of  the  model's  behavior.  The 
analytical  solution  is  obtained  for  normally  incident  waves  impinging  on  a 
plane  sloping  beach  and  is  used  to  compare  the  accuracy  of  the  numerical 
model.  The  analytical  solution  is  solved  by  integrating  (26)  for  the  wave 
height  (Thornton  and  Guza,  1983): 
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The  numerical  model  uses  a  modified  Euler  integration  scheme  to 
integrate  the  energy  flux  equation.  A  weighted  predictor/  corrector  loop  is 
used  to  increase  accuracy  and  stability  as  the  water  depth  becomes  shallow. 
The  basic  calculation  is  a  forward  stepping  scheme  with  a  variable  grid  size: 
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Well  outside  the  surf  zone  the  wave  height  changes  are  very  slight  with 
distance  as  the  bore  dissipation  is  essentially  zero,  and  are  relatively 
insensitive  to  the  grid  size.  However,  once  the  waves  start  to  break,  the 
bore  dissipation  increases  rapidly,  modifying  the  wave  heights.  To  increase 


the  resolution,  the  mesh  size  decreases  by  a  factor  of  four  when  the  bore 
dissipation  reaches  a  value  of  0.01.  This  flag  value  was  arrived  at  by  a 
series  of  model  tests  under  different  initial  conditions. 

The  model  was  tested  against  the  analytical  solution  for  a  variety  of 

initial  values.  The  initial  conditions  insured  that  shallow  water  criteria  were 

1  /2 

met  so  that  the  shallow  water  approximation  of  Cg  =  (gh)  ,  rather  than 
the  explicit  Equation  (24),  could  be  used.  The  first  test  used  a  shallow 
steepness  offshore  wave  (Hq/L0  =  0.01)  having  a  height  of  1.0  meter 
normally  impinging  on  a  plane  beach  with  a  slope  of  0.04.  The  period  of  the 
wave  was  10  seconds.  The  numerical  and  the  analytical  solutions  are 
compared  in  Figure  3.24.  Agreement  is  very  close  through  most  of  the 
offshore  grid  points.  As  the  water  depth  gets  very  shallow,  the  step  size  is 
again  reduced  to  insure  the  solution  converges  for  the  last  few  grid  points. 

For  a  steep  wave  case  (H  /L  =  0.02),  an  incident  wave  of  2.5  meters 

o  o 

was  used  as  the  offshore  input  (Fig.  3.24).  The  wave  period,  beach  slope 
and  normal  incidence  all  remained  the  same .  The  numeric  and  analytic 
solutions  are  again  in  good  agreement,  although  the  numeric  solution  slightly 
underpredicts  the  wave  height  as  it  did  in  the  case  where  Hq  equals  1.0  m. 
It  is  of  interest  to  note  that  the  high  wave  steepness  case  shows  no  growth 
in  wave  height  due  to  shoaling  prior  to  breaking.  The  wave  heights  decay 
monotonically  toward  the  shoreline. 

Having  satisfactorily  compared  the  model  results  with  the  analytic 
solution,  further  improvements  to  the  model's  computational  efficiency  were 
made.  Since  the  model  is  one  dimensional,  the  depth  variations  in  the 
longshore  direction  are  averaged  out.  Sand  bars  and  other  longshore  features 
are  taken  into  account  but  offshore  trenches  and  channels  are  smoothed  out 
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Figure  3.24.  Numeric  (solid)  versus  analytic  (dashed)  wave  height 
solutions  for  low  (II  =1.0  meter)  and  steep  (II  =2.5 
meter)  incident  wave  cases. 
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by  this  process.  The  new  height  module,  with  this  simplification,  runs  six 
times  faster  than  the  original  module. 

The  new  wave  height  subroutine  has  many  advantages.  It  uses  a  simple 
but  robust  probabilistic  approach  to  the  wave  height  problem.  It  allows 
variable  bathymetry  in  the  offshore  direction  and  obliquely  impinging  wave 
fronts.  It  is  fast,  accurate,  and  the  governing  equations  can  easily  be  fine 
tuned  by*  varying  the  parameters. 

D.  LONGSHORE  CURRENT  CALCULATION 

The  original  longshore  current  model  (Chen  and  Wang,  1983)  used  a  two 
dimensional  momentum  flux  balance  which  was  numerically  intensive.  The 
model  converged  very  slowly  (over  23  hours  of  CPU  time  in  some  cases),  and 
the  calculations  exhibited  unstable  behavior  (Devendorf,  1985).  To  provide 
the  necessary  longshore  current  information  to  an  operational  planner,  a 
simple  one  dimensional  model  was  substituted  for  the  original  two  dimensional 
subroutine.  The  model  was  designed  to  work  with  the  one  dimensional  wave 
height  model  discussed  in  the  previous  section.  The  random  wave  probability 
model  is  used  as  input  for  the  longshore  current  calculations.  A  narrow 
banded  wave  distribution  is  assumed  so  that  the  wave  frequencies  can  be 
approximated  by  a  mean  frequency  f  and  direction  i  .  The  longshore 
momentum  equation  is  a  balance  between  the  mean  wave  momentum  flux 
divergence  and  the  longshore  component  of  the  bed  shear  stress  (Thornton 
and  Guza,  1985): 
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b 
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(29) 


where : 


time  and  depth  averaged  covariance 
between  unsteady  velocity  components 

bed  shear  stress  in  longshore  direction 


S  is  made  up  of  two  terms  which  are  assumed  to  be  statistically 
yx 

*** 

independent.  The  first  term  is  the  radiation  stress  term  (S  )  which  is  the 

y  ^ 

wave  induced  momentum  flux.  The  second  term  (S  ')  is  the  depth 

yx 

integrated  turbulent  Reynold's  stress,  often  parameterized  in  terms  of  an 
eddy  viscosity  term. 

Monochromatic  current  models  require  the  eddy  viscosity  term  to  smooth 

out  the  discontinuity  in  the  current  field  at  the  breaker  line.  Since  no 

waves  break  outside  the  surf  zone  in  a  monochromatic  model,  there  is  no 

current  generated.  However,  there  is  a  jump  discontinuity  at  the  breaker 

line  which  requires  smoothing  by  the  eddy  viscosity  term.  In  the  random 

wave  model,  there  is  no  clearly  defined  breaker  line.  Some  waves  break 

offshore  and,  as  the  beach  is  approached,  more  waves  break  until,  in  the 

inner  surf  zone,  almost  all  waves  are  breaking  (Thornton  and  Guza,  1985). 

The  smooth  transition  of  energy  dissipation  reduces  the  requirement  for  the 

eddy  viscosity  term,  which  simplifies  the  velocity  calculation. 

The  radiation  shear  stress  S  is  described  by  (Longuet-Higgins  and 

yx 

Stewart,  1964): 


yx 
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Snell's  law  requires  that: 


S1-^-~  =  constant  =  - ^  (31) 

o 

so  that  when  (30)  is  differentiated  the  result  is: 

d  ;  sin  d  _  sin  eb 

dX,  Syx  ~  C  dX(E  Cg  cos  9)  =  C  <Cb>  (32) 

o  o 


This  implies  that  the  shear  stress  is  balanced  by  the  bore  dissipation  for 
breaking  waves  on  a  straight  and  parallel  beach. 

The  bed  shear  stress  is  given  as: 
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b 
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(33) 


where: 

[u|  =  instantaneous  total  bottom  velocity 
V  =  longshore  current  velocity 

The  bed  stress  is  linearized  by  assuming  small  angles  of  wave  incidence  and 
weak  longshore  currents  (Thornton,  1970;  Longuet-Higgins,  1970): 
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Using  shallow  water  wave  theory,  the  bottom  velocity  is  computed  by: 
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(35) 


Making  the  appropriate  substitutions  into  the  longshore  momentum  relation  and 
solving  for  V  s 

V 

The  numerical  scheme  is  compared  to  the  analytical  solution,  which  is 

obtained  -by  simply  substituting  the  expression  for  Hrmg  (Equation  27)  into 

the  velocity  calculation  (Equation  36).  The  numerical  and  the  analytical 

solutions  for  the  shallow  and  steep  wave  cases  studied  above  are  compared  in 

Figures  3.25  and  3.26.  The  numerical  solution  of  V  slightly  underpredicts 

the  current  velocity  compared  with  the  analytic  solution.  This  is  due  in 

part  to  the  introduction  of  a  small  error  in  H  which  was  discussed  in  the 

rms 

previous  section.  The  value  of  Hrmg  is  raised  to  the  sixth  power  in  Equation 
(36)  which  means  that  V  may  be  sensitive  to  small  errors  in  H 

rms 

However,  the  magnitude  of  the  error  for  V  is  on  the  order  of  0.05  and  0.07 
m/sec  which  may  be  neglected. 

The  initial  runs  of  the  updated  wave  height  and  longshore  current 
subroutines  are  presented  and  the  reworked  SURFCON  module  is  discussed  in 
the  next  section. 

E.  IMPROVEMENTS  TO  THE  SURFCON  MODULE 

During  the  reprogramming  effort,  a  major  consideration  was  maintaining 
the  integrity  of  the  original  model,  while  incorporating  modular 
improvements.  The  new  height,  direction  and  current  subroutines  were  coded 
under  new  names,  leaving  the  old  subroutines  intact.  To  maintain 
consistency,  the  control  module,  SURFCON,  was  renamed  OLDSURFCON  and 
a  new  module,  NEWSURFCON,  was  written  to  make  use  of  the  new 
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Figure  3.25.  Numeric  (solid)  versus  analytical  (dashed)  lonshore 
current  solutions  for  11  =1.0  meter.  Numerical  values 

are  slightly  lower  because  the  wave  height  values  are 
raised  to  the  seventh  power  in  the  velocity  calculation. 
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Figure  3.26.  Numeric  (solid)  versus  analytical  (dashed)  lonshore 
current  solutions  for  H  =2.5  meter. 


subroutines.  The  advantage  of  this  approach  is  that  the  model  can  be 
changed  from  the  new  to  the  old  version  simply  by  changing  the  main 
program  "CALL  NEWSURFCON"  to  "CALL  OLDSURFCON". 

The  WEWSURFCON  module  is  basically  the  same  as  the  old  module  with 
some  exceptions.  The  wavenumber  k,  celerity  C  and  group  speed  Cg  are 
initially  calculated.  Then  the  new  wave  direction  and  height  subroutines  are 
called,  which  return  direction,  height  and  longshore  current  matrices. 

The  first  major  change  incorporated  in  this  new  module  is  the  way  the 
breaker  line  is  defined.  In  the  old  module,  a  quadratic  equation  was  solved 
for  the  location  of  the  first  breaker  line.  Since  the  new  breaker  heights  are 
calculated  using  a  weighted  probability  density  function,  pb(H)  =  W(H)  p(H), 
there  is,  by  definition,  no  specific  spot  where  all  waves  are  said  to  be 
breaking.  The  weighting  function,  W(H),  as  defined  in  Equation  (37),  gives 
a  measure  of  the  percent  of  breaking  waves  (Thornton  and  Guza,  1983): 


W  (H) 


/  pb(H)d(H) 


(37) 


To  specify  a  "breaker  line" ,  the  model  first  differentiates  between  steep 
waves,  which  show  no  wave  growth  due  to  shoaling  towards  the  shore,  and 
shallow  waves,  which  increase  in  wave  height  until  they  break.  The  critical 
wave  steepness,  determined  from  Equation  (27),  is  Hq/Lo  =  0.02.  A  wave 
with  a  steepness  greater  than  this  criteria  is  considered  steep  while  a  wave 
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breaker  height".  That  is,  once  33  percent  of  the  waves  are  breaking,  the 
model  chooses  that  point  along  the  grid  as  the  first  breaker  line.  This  test 
is  actually  conducted  in  the  wave  height  subroutine.  As  the  model  moves 
down  the  column,  the  parameter  W(H)  is  tested  at  each  grid  point  to  see  if 
it  greater  than  or  equal  to  0.33.  Once  this  value  is  reached,  the  breaker 
height,  breaker  depth,  and  distance  from  shore  are  returned  to  the 
NEWSURPCON  module.  For  shallow  waves,  the  model  selects  the  maximum 
wave  height  as  the  breaker  line.  These  two  methods  give  good  estimates  of 
the  observed  breaker  line,  as  discussed  in  Chapter  5. 

Once  the  initial  breaker  height  is  defined  in  this  manner,  the  number  of 
additional  lines  are  computed  by  dividing  the  initial  breaker  distance  by  the 
wave  length  corresponding  to  the  average  surf  zone  depth.  If  there  are  one 
or  more  wave  lengths  between  the  initial  breaker  line  and  the  beach,  the 
number  of  breaker  lines  are  incremented  accordingly.  No  additional 
calculations  are  performed  on  these  interior  breaker  lines,  if  present.  The 
breaker  type  is  calculated  based  on  the  value  of  the  surf  parameter,  Eta, 
which  is  a  function  of  the  beach  slope  and  the  deep  water  wave  steepness 
(Battjes,  1974): 
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If  Eta  is  less  than  or  equal  to  0.4,  the  breaker  type  is  defined  as  spilling. 
Eta  values  between  0.4  and  2.0  are  defined  as  plunging  and  values  greater 
than  2.0  are  called  collapsing  or  surging  waves. 


The  next  calculation  performed  by  NEWSURFCON  is  the  "Effective  surf" 
computation.  Effective  surf  is  a  measure  of  surf  intensity  and  is  used  by 
operational  planners  to  provide  a  criterion  for  the  feasibility  of  conducting 
amphibious  landing  operations  using  specific  kinds  of  equipment.  Effective 
surf,  expressed  as  a  wave  height  in  feet,  is  a  modification  of  significant 
breaker  height,  taking  into  account  the  total  beach  and  surf  conditions.  It  is 
important  £0  note  that  effective  surf  is  simply  a  planning  parameter,  and  not 
a  direct  correspondence  to  the  significant  breaker  height.  Most  commonly, 
Hpff  is  larger  than  Hg,  but  in  some  cases,  especially  when  the  Hg  is  small 
and  the  breaking  waves  are  all  spilling,  Hpff  can  be  smaller  than  Hg. 

The  calculation  of  effective  surf  uses  a  simple  "look-up  table"  method 
which  is  outlined  in  the  Joint  Surf  Manual  ( COMNAVSURFPAC/ 
COMNA VSURFLANT  Inst.  3840.1,  1976).  The  calculations  of  effective  surf 
are  retained  from  the  original  model  with  only  minor  modifications. 
Originally,  the  effective  surf  was  retained  as  a  wave  height  matrix. 
However,  is  simply  a  single  number  for  a  specific  beach,  based  on 

significant  wave  height  and  direction,  period,  wind  speed  and  direction,  and 
percentage  of  spilling  waves.  The  new  module  calculates  a  single  number  for 
Hpff  and  outputs  all  the  parameters  as  specified  by  the  Joint  Surf 

Manual.  Other  information  of  use  to  the  planner  is  printed  at  the  same 
time.  The  new  calculation  also  allows  for  a  wind  input  (speed  and 
direction)  which  was  not  available  previously. 

The  final  operation  in  NEWSURFCON  is  the  limitation  of  the  wave  heights 
to  the  maximum  stable  wave  height  for  a  given  depth  and  period.  This 
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maximum  wave  height  is  given  as  (Wang  and  Chen,  1983): 


Hmax  =  0- 067  tanh  (kh)  (^)  (39) 

Each  generated  wave  height  is  compared  against  this  maximum  value  and  is 
reset  to  if  it  exceeds  the  limiting  wave  height.  The  wave  height, 

direction,  current  and  effective  surf  matrices  are  then  sent  to  OUTSURF  for 
formatted  output  to  either  the  screen  or  to  a  designated  printer. 

F.  DEVELOPMENT  OF  AN  OPERATIONAL  USERS  GUIDE 

A  simple  users  guide  was  written  for  the  new  version  of  the  SSSP  (Gill, 
1985).  Several  guidelines  were  adhered  to  during  its  development.  The 
manual  assumes  some  familiarity  with  the  HP-9845/275  and  a  working 
knowledge  of  oceanographic  terms.  The  manual  is  tutorial  in  nature.  It  is 
divided  into  two  parts,  for  the  two  main  modules  OPEN  and  SURFCON.  The 
manual  explains  each  input  request  and  appropriate  examples  are  included 
where  necessary.  An  effort  was  made  to  keep  the  manual  brief  and  user 
friendly.  The  model,  being  extensively  menu  driven,  is  straight-forward 
enough  to  require  little  documentation,  after  it  has  been  used  once  or  twice. 
The  BASIC  code  is  extensively  documented  throughout.  Where  appropriate, 
the  old  program  lines  are  left  in  place  (commented  to  prevent  execution), 
with  the  replacement  lines  documented  to  point  out  the  reason  for  the 
specific  change . 
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IT.  DISCUSSION 


A.  SENSITIVITY  TEST  OF  BREAKER  PARAMETERS 

As  discussed  in  Chapter  III,  the  computation  of  Hrmg,  using  the  Rayleigh 
probability  density  function  model,  depends  on  two  breaker  parameters, 
Gamma  and  B.  Gamma  is  an  adjustable  coefficient  of  0(1)  which  appears  to 
be  strongly  related  to  the  beach  slope  (Sallenger  and  Holman,  1985).  B  is 
the  breaker  coefficient  which  represents  the  percentage  of  foam  on  the  face 
of  the  breaking  wave  (i.e.,  the  intensity  of  the  breaking  wave).  This 
parameter  is  expected  to  be  less  than  or  equal  to  1.0,  since  a  B  value  of 
1.0  corresponds  to  a  fully  developed  bore  (a  wave  that  has  fully  broken  from 
crest  to  trough).  Values  of  B  greater  than  1.0  imply  that  the  bore 
dissipation  function  underestimates  the  wave  dissipation  due  to  breaking 
(Thornton  and  Guza,  1983). 

To  verify  the  validity  of  the  model  assumptions,  24  wave  data  sets  were 

studied  to  gain  an  understanding  of  the  interaction  of  the  B  and  Gamma 

values  as  used  in  the  model.  These  data  sets  are  summarized  in  Table  4.1. 

The  data  were  digitized  and  appropriate  Gamma  values  chosen,  based  on  the 

slope  of  H  /h  inside  the  surf  zone.  The  Gamma  values  ranged  from  0.39 
rms 

to  0.60  and,  for  each  data  set,  the  specific  Gamma  value  was  held  constant 
while  the  B  values  were  varied.  A  curve  fitting  routine  was  used  to 
optimize  the  B  parameter  for  a  given  Gamma.  As  B  was  varied,  the  new 
model  curve  was  compared  with  the  real  data  points,  and  a  sum  of  the 
square  of  the  errors  was  generated.  The  B  values  were  iterated  until  the 
error  sum  was  minimized.  Plots  of  the  model  output  versus  the  data  points 


were  generated  for  each  case. 


A  total  of  24  data  sets  were  tested  against  the  model  (See  Table  4.1). 

In  every  case,  the  model  accurately  predicts  the  wave  heights  described  by 

the  data  points.  The  first  ten  data  sets  are  from  a  laboratory  study  of 

energy  saturation  of  irregular  waves  during  shoaling  (Vincent,  in  press). 

These  wave  tank  tests  were  conducted  with  a  plane  bottom  slope  of  1:30, 

with  various  frequency  and  incident  wave  heights.  The  experiments  can  be 

divided  up  into  two  groups  based  on  wave  steepness.  The  low  steepness 

cases  (Hq/Lo  <  0.02)  are  represented  by  Case  1638  (Fig.  4.1).  The 

Gamma  value  of  0.60  was  derived  from  the  slope  of  H  /h  as  described 

rms 

above.  Using  the  curve  fitting  routine,  an  optimum  B  value  of  0.82  was 

chosen.  Note  the  good  fit  of  the  model  curve  (solid)  compared  with  the 

data  points  (circles).  The  high  steepness  cases  are  represented  by  Case 
1148  (Fig.  4.2).  Gamma  was  set  equal  to  0.60  and  a  B  value  of  1.2  was 
obtained.  As  before,  there  is  good  agreement  between  the  model  curve  and 
the  data  points. 

To  determine  the  optimized  model's  behavior  when  used  with  more 
difficult  bathymetry,  data  from  several  wave  tank  experiments  with  offshore 
bars  were  tested  (Battjes  and  Stive,  in  press).  Case  5  (Fig.  4.3)  has  a 
simple  offshore  bar  with  a  mean  bottom  slope  of  1:25.  Gamma  was  estimated 
to  be  .45  and  a  B  value  of  1.2  was  obtained.  Although  the  presence  of  the 
offshore  bar  (Fig.  4.4)  radically  changes  the  wave  height  data  profile,  the 
optimized  model  curve  again  faithfully  predicts  the  wave  heights.  Wave  tank 
test  case  15  (Fig.  4.5)  has  a  very  complex  bathymetry  (Fig.  4.6).  The 
model  divides  the  on /offshore  distance  into  10  grid  points  with  the 
intermediate  depths  calculated  by  straight  linear  interpolation.  Therefore, 
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DLE  4.1.  SURF  PAR 
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Figure  4.1.  Optimized  model  output  for  Case  1638  (Vincent,  1985). 


Optimized  model  output  for  Case  1148  (Vincent,  1985) 


OFFSHORE  DISTANCE  (CM) 

Optimized  model  output  for  Case  5  (Battjes  and  Stive, 


Optimized  model  output  for  Case  15  (Battjes  and  Stive, 


Figure  4.6.  Case  15  wave  tank  bathymetry.  (Battjes  and  Stive,  1985). 


complex  bathymetry,  such  as  that  in  Case  15,  was  smoothed  out  to  a 
certain  degree.  Even  with  that  limitation  in  mind,  the  optimized  model 
curve  appears  to  accurately  estimate  the  recorded  wave  heights.  The  model 
results  will  be  tested  against  field  data  in  the  next  chapter. 

It  is  evident  that  the  random  wave  model,  with  the  appropriate  choice 
of  the  surf  parameters,  Gamma  and  B,  will  accurately  forecast  wave  heights 
throughout  the  domain.  A  sensitivity  test  for  the  choice  of  parameters  was 
conducted  on  three  representative  cases  1638,  1148,  and  5.  Case  5  (Battjes 
and  Stive,  in  press)  was  chosen  for  discussion  here  because  of  its  more 
interesting  bottom  profile.  It  was  felt  that  this  barred  profile  was 
representative  of  many  real  beaches  and  would  provide  a  more  rigorous  test 
than  a  simple  plane  bottom.  Using  a  Gamma  value  of  0.45,  the  optimum  B 
value  of  1.3  was  changed  to  1.7,  1.5,  1.1  and  0.9  (Fig.  4.7).  B  values 
less  than  optimal  tend  to  over-estimate  the  wave  heights  and  vice-versa.  It 
is  of  interest  to  note  that  changes  in  B  of  22-33  percent  only  resu’ts  in  an 
increase  in  the  model  error  of  about  10-12%.  This  was  confirmed  in  a  study 
of  Torrey  Pines,  California  wave  data  by  Thornton  and  Guza  (1983).  They 
found  that  a  variation  of  +/-  25%  about  the  optimum  B  value  resulted  in  an 
increased  model  error  of  less  than  10%  (Fig.  4.8). 

It  was  found  that  the  model  demonstrated  a  similar  sensitivity  to  changes 
in  the  Gamma  value.  A  change  of  about  +/-  20%  to  the  optimum  Gamma 
value  of  0.45  resulted  in  an  increased  model  error-  of  8-10%  (Fig.  4.9). 
Sensitivity  tests  were  conducted  on  high  and  low  steepness  plane  beach  cases 
(1638  and  1148)  with  similar  results.  In  the  next  section,  techniques  and 
rationale  for  choosing  B  and  Gamma  values,  based  on  initial  conditions,  are 
discussed.  It  is  important  to  bear  in  mind  the  above  sensitivity  analysis. 
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Figure  4.7.  Sensitivity  test;  G  =  0.45;  B  varies  from  0.9  to 
about  the  optimum  B  value  of  1.3. 


ensitivity  test;  B  =  1.3;  Gama  varies  from  0.35  to 
.55  about  the  optimum  Gama  value  of  0.45.  The  model 
ppears  to  be  slightly  more  sensitive  to  flucuations  in 


The  model  appears  to  be  fairly  robust  with  respect  to  small  changes  in  the 
surf  parameters.  Therefore,  slight  errors  in  the  choice  of  these  values 
should  not  degrade  the  final  performance  of  the  model  to  a  great  extent. 

B.  CHOOSING  THE  BREAKER  PARAMETERS 

The  two  breaker  parameters,  B  and  Gamma,  are  multiplicative  constants 
in  the  Hrmg  wave  height  calculation  (see  Equation  27).  They  could 
therefore,,  be  combined  as  a  single  bulk  modulus  coefficient,  or  kept  as  two 
separate  parameters.  It  was  felt  that  keeping  the  two  terms  separate  would 
better  serve  the  requirements  of  the  model. 

Some  work  has  already  been  done  on  the  problem  of  choosing  the  breaker 
parameters.  A  recent  study  by  Sallenger  and  Holman  (1985)  has  determined 
a  functional  relationship  between  Gamma  and  the  average  beach  slope  based 
solely  on  field  observations  inside  the  surf  zone.  This  relationship  is  given 
by: 


Y  =  3.2  tan  8  +  0.3 


(40) 


Using  the  24  wave  data  sets,  the  chosen  Gamma  values  are  plotted  against 
beach  slope  (Fig.  4.10).  The  solid  line  is  a  straight  linear  regression  of  the 
data  points  while  the  dashed  line  is  the  Holman /Sallenger  function.  There 
seems  to  be  a  general  agreement  in  the  slope  of  the  line  although  the 
intercept  is  slightly  different.  Gamma  is  also  plotted  against  wave  steepness 
(Fig.  4.11)  but  the  plot  is  less  linearly  correlated.  As  an  initial  trial,  the 
Holman /Sallenger  functional  relationship  will  be  adopted  as  the  model's 
convention  for  choosing  the  Gamma  value.  The  average  slope  is  defined  here 
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as  the  depth  at  the  outermost  grid  point  divided  by  the  offshore  grid 


distance . 

Choosing  a  proper  predictor  for  the  B  value  proved  to  be  more  difficult. 
Since  B  represents  a  measure  of  breaker  intensity ,  it  is  felt  that  it  would  be 
related  to  incident  wave  steepness  which  has  previously  been  shown  to  be 
related  to  breaker  type  (Battjes,  1972).  The  optimum  B  values  for  the  24 
wave  sefs  are  plotted  against  steepness  in  Figure  4.12  and  against  beach 
slope  in  Figure  4.13.  The  wave  steepness  in  this  case  was  defined  as  the 
incident  wave  height  divided  by  the  deep  water  wave  length  (Ho/Lq). 
Although  there  is  a  considerable  scatter,  B  appears  correlated  with  wave 
steepness  and  less  with  beach  slope. 

Several  other  parameters  were  considered  as  predictors  for  B.  The  B 

and  Gamma  values  were  compared  with  the  surf  parameter,  Eta  =  tan  8  / 

3  4 

(Ho/Lq)  .  Additionally,  the  Bulk  Modulus  (B  /Gamma  )  was  plotted 
against  beach  slope  and  wave  steepness.  Since,  from  Figures  4.11  and  4.12, 
B  appears  to  be  related  to  both  beach  slope  and  wave  steepness,  it  was 
hoped  that  Eta,  which  is  a  measure  of  the  expected  breaker  intensity,  would 
offer  a  predictor  for  B.  The  plot  of  B  versus  Eta  (Fig.  4.14)  shows  an 
interesting  pattern.  B  appears  to  increase  almost  linearly  as  Eta  increases 
to  about  0.6  and  then  decreases  in  an  exponential  fashion  to  an  asymptotic  B 
value  of  0.8.  This  behavior  appears  reasonable.  As  the  surf  parameter 
increases,  indicating  a  transition  from  spilling  to  plunging  waves,  the  B  value 
( amount  of  wave  face  that  is  breaking ) ,  should  increase . 

The  maximum  value  that  B  should  reach,  according  to  theory,  is  1.0. 
Studies  have  shown,  however,  that  in  many  cases,  the  optimum  B  value  is 
on  the  order  of  1.4  (Thornton  and  Guza,  1983).  The  maximum  B  value 
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SURF  PARAMETER  (ETA) 

Figure  4.14.  B  versus  the  surf  parameter  Eta.  B  is  strongly  related  to 
Eta  which  is  to  be  expected.  Both  B  and  Eta  are  measures  of 
the  intensity  of  the  wave  breaking.  Dashed  line  depicts 
piecewise  function  used  by  the  model  to  define  B. 


occurs  when  Eta  is  about  0.6,  indicating  that  the  breakers  are  mixed  spilling 
and  plunging.  As  Eta  continues  to  increase,  the  B  values  fall  more  into  line 
with  theory  and  approach  the  expected  value  of  1.0,  as  the  waves  become 
fully  plunging. 

To  simplify  and  expedite  the  model  development,  the  following  piecewise 
continuous  function  is  used  to  determine  the  B  term  (Fig.  4.14): 

B  =  2.167  Eta  +  0.4  Eta  <=  0.6  (41) 

B  =  -  1.70  Eta  +  2.7  0.6  <  Eta  <=  1.0 

B  =  1.0  Eta  >  1.0 

This  choice  is  somewhat  ad  hoc  and  is  based  primarily  on  the  data  sets 
discussed  above.  Future  refinements  in  the  techniques  for  predicting  the 
breaker  parameters  may  substantially  improve  the  use  of  this  model 
operationally.  Using  equations  (40)  and  (41)  as  predictors  for  Gamma  and 

B,  the  SURFCON  module  is  essentially  complete.  The  next  chapter  will 
discuss  testing  the  model  against  two  well-documented  natural  beach  data 
sets,  to  get  a  feeling  for  how  the  improved  model  behaves  against  actual 
data. 

C.  MOTIVATION  FOR  THE  RANDOM  WAVE  MODEL  SUBSTITUTION 

As  discussed  in  Chapters  3  and  4,  the  primary  motivation  for  substituting 
the  one  dimensional  random  wave  model  for  the  two  dimensional  model  was  to 
increase  computational  efficiency  and  accuracy.  The  improved  model 
calculates  longshore  current  information  quickly  and  accurately.  Significant 
wave  height  predictions  compare  well  with  actual  data  sets,  as  will  be 
discussed  in  Chapter  5.  There  is  ,  however,  a  potential  for  loss  of  detailed 


information  by  using  the  one  dimensional  model.  Nearshore  circulation  and 
longshore  variation  in  wave  heights  are  smoothed  out  by  the  one  dimensional 
approach . 

The  surf  model  (both  old  and  new  versions)  uses  a  single  wave  height 
input  to  define  the  offshore  boundary  condition.  This  implies  that  variations 
in  the  interior  wave  heights  and  directions  are  due  only  to  variations  in 
bathymetry.  Because  the  model  uses  a  very  small  grid  (only  13  longshore 
grid  points),  very  accurate  bathymetry  must  be  entered.  It  is  questionable 
whether  the  operational  planner  would  have  access  to  bathymetric  data  of 
the  required  accuracy  to  define  a  two  dimensional  bathymetry. 

Many  beaches  which  would  potentially  be  used  for  amphibious  operations 
are  shallow,  almost  planar  in  nature.  Other  common  landing  beaches  are 
barred,  but  rather  uniform  in  the  longshore  direction.  Therefore,  for  many 
operationally  relevant  beaches ,  the  assumption  of  straight  and  parallel 
contours,  and  the  use  of  a  one  dimensional  model,  appears  justified.  The 
"mean  beach"  surf  information  provided  by  the  model  will  adequately 
represent  the  conditions  over  a  standard  500  yard  landing  beach.  If  a  wider 
beach  needs  to  be  characterized,  the  improved  model  would  be  run  several 
times,  entering  the  required  information  for  each  beach  segment. 
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TEST  AND  EVALUATION  OF  THE  IMPROVED  MODEL 


A.  WAVE  HEIGHT  AND  BREAKER  LOCATION 

The  improved  version  of  the  SSSP  was  tested  against  one  wave  tank 
experiment  and  two  sets  of  beach  surf  data  to  see  if  the  model  would  choose 
reasonable.  Gamma  and  B  parameter  values,  based  on  the  predictors  discussed 
in  Chapter  4.  The  predicted  wave  heights  and  surf  characteristics  are 
compared  with  the  observed  values. 

Case  5  (Battjes  and  Stive,  1985),  discussed  in  the  preceding  chapter, 
was  chosen  as  the  wave  tank  experiment  with  which  to  test  the  finished 
model.  This  case  was  chosen  because  of  its  barred  beach  bathymetry. 
The  model  selected  values  of  0.41  and  0.963  for  Gamma  and  B.  The 
optimum  values,  from  Table  4.1,  are  0.45  and  1.30.  The  model  wave  height 
fit  to  the  observed  data  points  is  still  quite  reasonable  (Fig.  5.1).  The 
solid  line  represents  the  model  output  and  the  circles  are  the  observed  data. 

Two  Southern  California  beach  data  sets  were  used  for  further  tests. 
Data  collected  at  Santa  Barbara's  Leadbetter  Beach  during  February,  1980 
and  Torrey  Pines  during  November,  1978  (Thornton  and  Guza,  1983),  were 
used  as  the  primary  tests  of  the  performance  of  SSSP. 

The  first  Santa  Barbara  test  used  data  collected  on  3  February,  1980. 
The  average  beach  slope  was  0.047  and  the  incident  waves  were  0.55  meters 
(rms)  with  a  period  of  14.3  seconds.  The  observed  breaker  height  was 
measured  at  0.70  meters,  approximately  40  meters  from  the  beach  (dashed 
line).  The  wave  angle  was  7.8  degrees  from  normal.  The  SSSP  predicted  a 
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Figure  5.1.  Comparison  between  the  model  output  (solid)  and  the 
observed  wave  tank  data  (dashed).  C,  for  all  cases  is 
0.009.  f 


rms  breaker  height  of  0.68  meters  (solid  line),  40.0  meters  from  the  beach 
(Fig.  5.2).  The  model  under  predicted  the  observed  heights  less  than  4 
percent . 

Similar  results  are  obtained  for  Leadbetter  Beach  data  obtained  on  4  and 
5  February,  1980.  On  4  February,  the  incident  wave  height  was  0.56 
meters  with  a  period  of  14.3  seconds.  The  wave  angle  was  9  degrees  from 
normal.  The  observed  breaker  height  was  0.67  meters,  50  meters  from  the 
beach.  The  model  predicted  0.65  meters,  40  meters  from  the  shore  (Fig. 
5.3) . 

On  5  February,  the  incident  wave  was  0.45  meters  with  a  period  of 
12.8  seconds.  The  wave  angle  was  8.4  degrees  from  normal.  Observed 
breakers  were  0.60  meters,  about  36  meters  from  the  beach.  SSSP  predicted 
0.55  meter  breakers,  33  meters  from  the  shore  (Fig.  5.4). 

Torrey  Pines  beach  data  acquired  on  4  and  10  November,  1978  is 
compared  with  model  results  in  Figures  5.5  and  5.6.  The  beach  has  a  shallow 
slope  of  0.022  which  accounted  for  a  wide  surf  zone  of  about  160  meters. 
On  4  November,  the  incident  waves  were  about  0.35  meters  with  a  period  of 
14.3  seconds  and  normal  incidence.  The  observed  breaker  height  was  0.47 
meters  about  110  meters  from  the  shore.  SSSP  predicted  0.48  meter 
breakers  at  a  distance  of  106  meters  from  the  beach.  The  model  correctly 
forecast  the  wave  height  and  surf  zone  width  (Fig.  5.5). 

A  similar  result  is  seen  for  the  November  10  Torrey  Pines  Case.  The 
incident  wave  is  0.68  meters  with  a  period  of  15.9  seconds  and  normal 
incidence.  Observed  breakers  were  0.72  meters,  approximately  130  meters 
from  the  shoreline.  The  model  predicted  breaker  heights  of  0.75  meters,  140 
meters  from  the  beach  (Fig.  5.6).  As  for  the  November  4  case,  the  model 
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Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Santa  Barbara  on  3  Feb., 
(dashed) . 


orrsnoriR  distance  (m) 

Figure  5.3.  Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Santa  Barbara  on  A  Feb.,  1980 
(dashed) . 


Figure  5.4,  Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Santa  Barbara  on  5  Feb.  ,  1980 

(dashed) . 


Figure  5.5.  Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Torrey  Pines  on  4  Nov.,  1978 
(dashed) . 


Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Torrey  Pines  on  10  Nov.,  1978 
(dashed) . 


worked  quite  well  for  wave  height  prediction  and  the  location  of  the  breaker 
line . 

Based  on  the  tests  of  the  SSSP  against  various  data  sets,  it  appears  that 
the  new  version  is  an  improvement  over  the  original  model.  The  old  version 
consistently  underpredicted  the  wave  heights  and  was  very  inaccurate  at 
estimating  the  location  of  the  breaker  line  (Devendorf,  1985).  The  new 
version  appears  to  be  quite  good  at  predicting  the  breaker  height  and  gives 
better  estimates  of  the  breaker  line.  Table  5.1  is  a  comparison  between  the 
optimum  values  of  the  two  breaker  parameters  and  the  values  chosen  by  the 
model.  The  technique  for  choosing  Gamma  and  B  is  convenient,  although 
improvements  will  need  to  be  made  in  this  area. 

B.  LONGSHORE  CURRENT  DISTRIBUTION 

Longshore  currents  are  the  result  of  waves  impinging  at  an  oblique  angle 
to  the  shoreline.  All  the  wave  tank  tests,  as  well  as  the  Torrey  Pines 
cases,  had  normally  incident  waves.  The  Santa  Barbara  Leadbetter  Beach 
data  is  the  only  comprehensive  field  data  on  longshore  currents  acquired  to 
date.  There  exists  no  laboratory  data  for  longshore  currents  generated  by 
obliquely  incident,  random  waves.  The  same  three  data  sets  (February 
3,4,5)  were  chosen  as  study  cases  because  they  had  a  variety  of  incident 
wave  heights  and  wave  angles.  For  the  February  3  case  (Fig.  5.7),  the 
maximum  observed  current  (dashed)  was  0.48  m/sec  at  a  distance  of  22 
meters  from  the  shoreline.  The  bed  shear  stress  coefficient,  C^,  was  set 
equal  to  0.009  for  all  cases.  This  value  was  derived  experimentally  by 
Thornton  and  Guza  (1985).  The  SSSP  prediction  (solid)  for  longshore  current 
was  a  maximum  of  0.5  m/sec,  25  meters  from  the  beach  (solid  line). 
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For  the  February  4  case  (Fig.  5.8),  the  maximum  observed  current  was 
.48  m/sec,  25  meters  from  shore;  the  model  prediction  was  .44  m/sec  at 
20  meters  from  the  beach.  The  Feb  5  case  (Fig.  5.9)  shows  an  observed 
current  maximum  of  .36  m/sec  at  20  meters  from  the  beach.  The  model 
predicted  .36  m/sec  at  20  meters  from  the  shoreline. 

The  model  not  only  predicts  the  maximum  current  velocity  reasonably 
well,  but  the  current  distribution  within  the  surf  zone  shows  good  agreement 
with  the  data.  The  model  seems  to  underpredict  the  current  velocity 
seaward  of  the  maximum  current,  but,  in  general,  the  model  does  quite  a 
good  job  of  longshore  current  prediction. 

The  other  Leadbetter  Beach  cases  showed  similar  results.  In  all  cues, 
the  model  gives  reasonable  maximum  current  locations  and  speeds  as  well  as 
current  distributions  within  the  surf  zones .  These  findings  were  also 
confirmed  by  Thornton  and  Guza  (1985)  in  a  study  in  which  the  random  wave 
model's  wave  height  and  longshore  current  predictions  are  rigorously  tested 
against  the  Leadbetter  Beach  data. 
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Figure  5.7.  Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Santa  Barbara  on  3  Feb.,  1980 
(triangles) . 


Figure  5.8.  Comparison  between  the  model  output  (solid)  and  the 
observed  data  from  Santa  Barbara  on  4  Feb.,  1980 
( triangles) . 


VI.  SUMMARY  AND  CONCLUSIONS 

A.  SUMMARY 

In  1981,  the  U.S.  Navy  contracted  to  have  a  Sea,  Swell  and  Surf 
Program  (SSSP)  developed  to  provide  input  to  the  planning  efforts  of  Naval 
Oceanographers  at  sea  and  amphibious  operational  planners.  The  model  was 
developed  in  FORTRAN  and  designed  to  run  on  a  main  frame  computer  (Wang 
and  Chen,  1983).  The  completed  SSSP  was  then  translated  into  BASIC  for 
use  on  the  currently  available  shipboard  micro-computers,  the  HP-9845B/275 
(Devendorf,  1985).  The  SSSP  was  modular  in  nature,  with  separate 
subroutines  for  offshore  wind  wave  generation ,  surf  zone  calculations  and 
nearshore  circulation. 

The  original  SSSP  had  several  problems  which  prevented  its  use  by 
operational  planners.  The  offshore  wind  wave  generation  module  overbuilt 
the  wave  heights  for  high  wind  speeds.  Some  of  the  surf  condition 
calculations  were  inaccurate  and  the  nearshore  circulation  subroutine  had 
stability  problems.  These  problems  were  addressed  by  this  thesis  in  an 
attempt  to  improve  and  upgrade  the  SSSP. 

The  open  water  module  shows  satisfactory  wind  rave  generation  for  wind 
speeds  of  15  knots  or  less.  For  higher  wind  speeds,  the  predicted  wave 
generation  occurs  too  quickly,  both  spatially  and  temporally.  Sensitivity 
analysis  indicated  that  the  Phillips  linear  term  was  responsible  for  the  rapid 
development,  but  no  improvements  to  the  model  physics  were  made.  Several 
ad  hoc  "Axes"  were  recommended.  This  problem  will  require  further 
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Two  subroutines  are  called  to  calculate  the  wave  transformation  as  the 
wave  proceeds  from  offshore  to  onshore.  Several  numeric  improvements  were 
made  which  eliminated  minor  instabilities  in  the  wave  direction  and  height 
subroutines.  The  original  no  flux  boundary  conditions  were  changed  to  cyclic 
boundary  conditions  in  both  subroutines.  This  change  allows  the  wave 
heights  and  directions  to  more  accurately  reflect  the  conditions  required  by 
the  model  physics  at  the  model  boundaries. 

The  most  serious  problem,  the  instability  of  the  nearshore  circulation 
module,  is  solved  by  substituting  a  one  dimensional  longshore  current  model 
for  the  original  two  dimensional  scheme.  Although  there  is  a  loss  of  two 
dimensional  current  resolution  within  the  surf  zone  due  to  the  assumptions  of 
straight  and  parallel  contours,  stationary  wave  conditions,  and  shallow  water 
initial  waves,  it  was  felt  that  a  model  that  gives  the  important  longshore 
current  information  quickly  and  accurately  is  an  improvement  over  the 
original  unstable,  two-dimensional  model.  The  new  current  model  is  a 
simple,  longshore  momentum  flux  balance  equation.  It  retains  the  wave 
induced  momentum  flux  and  neglects  the  turbulent  Reynold's  stress.  The  new 
model  is  fast  and  accurate,  requiring  only  a  few  calculations  once  the  wave 
height  field  is  generated. 

To  take  advantage  of  state  of  the  art  breaker  height  modeling,  a  random 
wave  model  was  implemented  for  the  calculations  of  the  wave  heights  in  the 
nearshore  area.  The  random  wave  model  assumes  that  wave  heights,  in  the 
nearshore  region,  follow  a  Rayleigh  distribution  and  that,  even  after  the 
waves  break,  this  distribution  is  retained. As  a  future  project,  the  random 
wave  model  could  be  developed  Into  a  two  dimensional  model  by  applying  the 
proper  longshore  numerical  differencing. 


The  random  wave  model  assumes  that  wave  breaking  can  be  described  by 
a  periodic  bore.  Two  parameters  are  introduced  by  the  model,  which  must 
be  chosen  based  on  the  initial  conditions.  The  Gamma  term  is  part  of  the 
proportionality  coefficient  which  modifies  the  initial  Rayleigh  distribution,  and 
is  found  to  be  strongly  correlated  with  the  beach  slope.  The  breaker 
parameter,  B,  describes  the  amount  of  the  wave  face  that  is  breaking  and 
appears  to  be  related  to  the  surf  parameter  Eta.  The  criterion  for  choosing 
the  two  breaker  parameters  was  developed  empirically. 

The  wave  height  and  current  models  were  next  tested  against  wave  tank 
and  natural  beach  data  sets.  Although  the  model  did  not  necessarily  choose 
optimum  values  for  the  breaker  parameters,  the  model  is  fairly  robust,  and 
slight  deviations  in  the  parameter  values  do  not  degrade  the  wave  height 
calculations  significantly.  Good  agreement  was  found  between  the  observed 
wave  heights  and  the  model's  predicted  values.  The  longshore  current,  based 
on  the  model's  predicted  wave  height  field,  were  also  in  close  agreement 
with  observations. 

The  calculation  of  Effective  Surf,  an  amphibious  planning  guideline,  was 
improved.  The  algorithm  was  implemented  in  accordance  with  the  Joint  Surf 
Manual  (1976)  with  respect  to  its  input  requirements  and  its  formatted 
output. 

A  simple  users  manual,  tutorial  in  nature,  was  developed  for  this  version 
of  the  SSSP.  It  is  short,  and  assumes  some  knowledge  of  the  host  computer 
and  basic  oceanographic  terms.  It  is  written  within  the  guidelines  set  forth 
by  the  Navy  for  development  of  computer  models  for  the  Shipboard  Numerical 
Aid  Program  (Brown,  1984). 
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B.  CONCLUSIONS 


The  improved  Sea,  Swell  and  Surf  Program  (SSSP)  can  provide  much 
useful  information  to  the  Naval  Oceanographer  at  sea  or  amphibious  planner. 
It  has  several  problems  which  have  not  been  resolved,  so  its  use  is 
recommended  with  several  caveats.  The  offshore  wind  wave  generation 
module  is  felt  to  be  accurate  for  low  wind  speeds  only.  Useful  information 
on  offshore  wave  growth  can  be  obtained  for  in^ut  wind  speeds  of  15  knots  or 
less.  The  results  obtained  for  higher  wind  speed  inputs  must  be  examined 
carefully  before  use.  The  requirement  for  providing  an  initial  wave  energy 
spectrum  is  probably  impractical  for  shipboard  use.  The  option  of  providing  a 
single  initial  wave  height  and  period  will  be  of  more  value  to  the  operational 
user. 

The  new  random  wave  breaker  height  model  appears  to  work  well.  The 
breaker  heights  and  locations  as  well  as  the  longshore  current  field  are 
accurately  predicted  by  the  model.  Due  to  its  one-dimensional  nature  and 
the  assumption  of  straight  and  parallel  contours,  longshore  wave  height  and 
current  resolution  is  restricted.  However,  the  assumption  that  an  amphibious 
landing  beach  is  uniform  in  the  longshore  direction  is  a  reasonable  one  over  a 
wide  range  of  operational  areas.  The  small  grid  size  of  the  model  supports 
the  use  of  a  one  dimensional  model  to  compute  an  average  set  of  charac¬ 
teristics  for  a  "mean"  beach.  The  speed  and  stability  of  the  new  model 
make  it  a  viable  tool  for  operational  planners.  Recommendations  for  further 
improvement  of  the  SSSP  include: 

1.  Pursue  the  problem  of  the  overbuilding  of  the  offshore  wind  waves. 

2.  Development  of  better  predictors  for  the  breaker  parameters  Gamma 
and  B  by  including  a  wider  data  base.  As  was  seen  in  the 
optimized  tests,  correct  choice  of  these  parameters  will  provide  a 
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very  good  fit  to  observations,  even  with  very  complicated 
bathymetry. 

3.  Explore  the  feasibility  of  using  the  random  wave  model  two 
dimensionally.  This  would  offer  the  advantages  of  this  more  realistic 
approach  to  wave  height  modeling,  while  retaining  the  details  that 
two  dimensional  spatial  resolution  provides. 

The  improved  SSSP  offers  the  Naval  Oceanographer  at  sea  and  amphibious 
operational  planners  a  new  and  useful  tool  that,  when  used  in  conjunction 
with  surf  observations,  weather  information,  and  other  available  data,  can  be 
a  great  aid  in  the  safe  planning  and  execution  of  beach  landing  operations. 
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